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eddies.  The  second  scale  is  based  on  the  well  known  Kolmogorov  hypothesis 
that  dissipation  of  turbulent  kinetic  energy  occurs  primarily  at  small  eddies. 
The  turbulence  model  derived  based  on  the  concept  of  two  different  scales 
is  called  the  two-scale  turbulence  model.  The  existing  turbulence  model 
which  is  modelled  based  on  the  one-scale  concept  of  k  and  e  is  called  the 
one-scale  turbulence  model. 

The  two-scale  turbulence  model  is  then  applied  to  predict  turbulent 
free  shear  flows  and  recirculating  flows.  The  calculations  were  done  in 
three  parts.  The  first  test  case  was  nonbuoyant  free  shear  flows  which 
included  round  and  plane  jets  in  stagnant  and  moving  streams,  plane  wakes  and 
mixing  layer.  In  the  second  part,  the  model  was  tested  for  plane  and  round 
buoyant  jets  having  different  Froude  numbers.  Finally,  some  results  were 
obtained  for  recirculating  flows,  namely,  backward  facing  step  and  flow  past 
an  obstruction. 

It  is  shown  in  the  present  study  that  the  two-scale  turbulence  model 
performs  significantly  better  than  the  one-scale  turbulence  model  in  all 
the  cases  concerned.  The  prediction  capability  of  the  two-scale  turbulence 
model  is  shown  since  one  does  not  need  to  alter  or  modify  the  turbulence 
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ABSTRACT 


The  use  of  second  order  closure  turbulence  model  in  predicting 
turbulent  flows  is  known  to  be  more  successful  than  the  classical  mixing 
length  model.  However,  it  is  found  that  if  the  turbulence  constants  are 
not  altered  or  modified,  the  second  order  closure  turbulence  model  is 
unable  to  predict  satisfactorily  for  some  flows  such  as  round  jet  and 
wake  flows.  In  order  to  improve  the  predictability  of  the  second  order 
closure  model,  the  present  work  proposes  to  consider  two  turbulent 
scales  in  the  modelling  of  turbulent  flows.  One  of  these  scales  is  based 
on  using  the  turbulent  kinetic  energy,  k,  and  its  dissipation  rate,  e, 
to  characterize  the  large  energy  containing  eddies.  The  other  scale  is 
based  on  the  dissipation  rate,  s,  and  the  kinematic  viscosity,  v,  to 
characterize  the  small  energy  dissipating  eddies.  The  second  scale  is 
based  on  the  well  known  Kolmogorov  hypothesis  that  dissipation  of 
turbulent  kinetic  energy  occurs  primarily  at  small  eddies.  The 
turbulence  model  derived  based  on  the  concept  of  two  different  scales  is 
called  the  two-scale  turbulence  model.  The  existing  turbulence  model 
which  is  modelled  based  on  the  one-scale  concept  of  k  and  s  is  called 
the  one-scale  turbulence  model. 

The  two-scale  turbulence  model  is  then  applied  to  predict 
turbulent  free  shear  flows  and  recirculating  flows.  The  calculations 
were  done  in  three  parts.  The  first  test  case  was  nonbuoyant  free  shear 


iii 


flows  which  included  round  and  plane  jets  in  stagnant  and  moving 
streams,  plane  wakes  and  mixing  layer.  In  the  second  part,  the  model  was 
tested  for  plane  and  round  buoyant  jets  having  different  Froude  numbers . 
Finally,  some  results  were  obtained  for  recirculating  flows,  namely, 
backward  facing  step  and  flow  past  an  obstruction. 

It  is  shown  in  the  present  study  that  the  two-scale  turbulence 
model  performs  significantly  better  than  the  one-scale  turbulence  model 
in  all  the  cases  concerned.  The  prediction  capability  of  the  two-scale 
turbulence  model  is  Shown  since  one  does  not  need  to  alter  or  modify  the 
turbulence  constants  as  in  the  case  of  the  one-scale  turbulence  model. 
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axisymmetric  case 
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CHAPTER  I 
INTRODUCTION 

1.1  Motivation  of  Research 

Many  fluid  motions  that  occur  in  nature  are  turbulent,  e.g.  flow  over 
aeroplanes,  ships  and  cars,  flow  in  jet  engines  and  turbines,  flow 
through  pipes  and  ducts ,  weather  patterns  and  ocean  waves .  Turbulent 
flow  is  a  complex  phenomena  that  plays  an  important  role  in  many 
engineering  designs.  Therefore,  it  is  important  for  engineers  to  study 
and  understand  this  complex  flow  and  be  able  to  predict  it.  Equations 
for  describing  the  fluid  motions,  known  as  the  Navier-Stokes  equations, 
have  been  postulated  and  derived  for  over  a  century.  However,  it  is 
difficult  to  solve  these  equations  for  both  laminar  and  turbulent  flows 
mainly  due  to  the  nonlinearity  of  the  equations.  For  turbulent  flows, 
the  problem  is  even  more  formidable  because  the  turbulent  fluid  motion 
is  irregular ,  random,  time  dependent  and  three  dimensional.  However,  in 
most  engineering  applications,  the  detailed  analysis  of  instantaneous 
turbulent  motion  is  not  necessary  and  the  gross  parameters  like  mean 
velocity,  average  pressure,  mean  temperature,  wall  shear  stress  and  wall 
heat  flux  are  often  sufficient  for  engineering  design. 

In  1895,  0.  Reynolds  [1]  proposed  an  averaging  technique  by  assuming 

JL. 

that  the  variable  0  at  any  instant  consists  of  the  mean  quantity  0  and 
a  fluctuating  part  0'.  Hence, 
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0  —  0  4-  0  1 

The  time  averaging  process  when  applied  to  the  Navier-Stokes  equations, 
creates  six  additional  unknowns.  These  unknowns,  although  called 
Reynolds  stress,  are  created  from  the  convective  or  non-linear  terms  of 
the  Navier-Stokes  equations.  In  order  to  solve  the  turbulent  flow 
problem  from  the  time  averaged  Navier-Stokes  equations  more  equations  or 
empirical  relations  are  needed  for  Reynolds  stress.  Methods  for 
deriving  equations  which  specify  a  relation  between  the  Reynolds  stress 
and  the  mean  flow  quantities  are  called  turbulence  models.  In  other 
words,  a  turbulence  model  is  needed  to  recover  the  information  of 
turbulent  motion  that  is  lost  in  the  averaging  process.  There  are  many 
turbulence  models  proposed  to  date.  However,  these  models  can  predict 
accurately  time  averaged  turbulent  flows  only  for  a  certain  class  of 
problem.  A  more  general  model  is  needed  if  one  expects  a  turbulence 
model  to  have  a  better  prediction  capability  and  a  practical  value  for 
engineering  applications. 

The  purpose  of  this  research  work  is  to  introduce  a  new  physical 
concept  into  the  modelling  of  turbulent  flows  and  to  improve 
predictability  of  the  model.  The  new  model  is  developed  for  the  second 
order  turbulence  correlation  based  on  the  concept  of  two  turbulent 
scales,  one  for  large  or  energy  containing  eddies  and  the  other  for 
small  or  energy  dissipating  eddies.  The  two-scale  turbulence  model  is 
first  tested  and  verified  for  a  class  of  turbulent  flows  called  'Free 
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Turbulent  Shear  Flows'.  In  free  turbulent  shear  flows,  shear  stress, 
neat  flux  and  diffusion  are  significant  in  the  directions  perpendicular 
to  the  direction  of  flow  and  there  is  no  solid  wall  in  the  flow  domain. 
Some  examples,  shown  in  figure  1.1,  are  mixing  layer,  coaxial  and  plane 
jets,  plumes,  buoyant  jets  and  wakes.  The  two-scale  turbulence  model  is 
then  used  to  predict  some  turbulent  separation  phenomena  such  as  flow 
separation  behind  a  step  as  shown  in  figure  1.2. 

There  are  several  reasons  for  selecting  free  turbulent  shear  flows  to 
test  the  turbulence  model.  First,  free  shear  flows,  as  shown  in  figure 
1.1,  have  a  weak  pressure  gradient  so  that  the  flow  characteristic  is 
largely  controlled  by  turbulent  shear  motion  which  affects  diffusion, 
production  and  dissipation  of  turbulent  motion  and  not  by  pressure 
force.  Therefore,  the  prediction  of  turbulent  free  shear  flow  is  more 
sensitive  to  the  turbulence  model  than  flows  with  large  pressure 
gradient.  Secondly,  abundant  experimental  data  are  available  and 
comparison  between  predicted  and  experimental  results  can  be  made  in 
detail.  Thirdly,  the  complication  of  near  wall  turbulence  is  not 
present  in  free  shear  flows  so  that  the  accuracy  of  the  two-scale 
turbulence  model  can  be  carefully  examined  without  the  interference  of 
wall  turbulence.  Fourthly,  turbulent  shear  flows  have  a  number  of 
practical  applications  and  play  an  important  role  in  various  engineering 
design.  Jet  engines,  chimney  plumes,  jet  streams  in  atmosphere,  wakes 
behind  aeroplanes  and  ships  and  cooling  water  disposal  in  rivers  are 
some  of  the  examples.  Though  some  of  these  flows  have  walls  in  their 
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Figure  1.1.  Examples  of  Free  Shear  Flows 
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Figure  1.2.  Examples  of  Separated  Flow 


vicinity,  the  study  of  free  shear  flows  is,  nevertheless,  a  first  step 
in  understanding  problems  and  phenomena  involved. 
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1.2  Historical  Development  of  Turbulence  Models 


In  this  section,  a  brief  historical  review  of  turbulence  modelling  is 
made  leading  to  a  discussion  of  the  problems  in  some  of  the  models.  In 
order  to  resolve  the  difficulties  in  the  existing  models,  a  new  model  is 
presented. 

As  mentioned  earlier,  the  need  of  turbulence  modelling  arose  when 
Reynolds  [1]  proposed  the  averaging  process  to  obtain  governing 
equations  for  turbulent  flows.  To  faciltate  the  discussion,  the- Navier- 
Stokes  equations  and  the  energy  equation  for  incompressible  flow  are 
written  here  as 


(1.1) 


_ i 

Dt 


p3x 


3P 


(1.2) 


*  ic  * 

DT  *  3U.  32T 


(1.3) 


J  j  j 


The  instantaneous  quantities  for  velocity,  pressure,  stress  and 


*  -fr  iz  V? 

temperature  ^  ,  P  ,  ,  T  are  denoted  by 


U. 

1 
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where  the  quantities  on  the  right  are  the  mean,  U  ,  P,  t  T  and 

i  ij’ 

fluctuation,  u  ,  p,  T..',  0,  of  velocity,  pressure,  stress  and 

1  J 

temperature.  These  are  substituted  in  the  Navier-Stokes  equations  and 
averaged  by  a  short  time  average  or  ensemble  average  to  give 


* 

3U. 

_ i 

3x . 


=  0 


(1.4) 
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Dt 


3P 
p  3x . 
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3ZU . 


+  v 


3x . 3x . 
J  J 


3u  ,u . 
_ L  J 

3x . 

J 


(1.5) 


pc 


.DT 

dF 


3U.  3ZT  3u . 0 

-  Tii _ 1  +  zr _ i_ 

J3x.  3x.3x.  3x . 

J  J  J  i 


t . . 1 3u . 

+  -il  T-i 
pc  a*j 


(1.6) 


These  set  of  equations  introduce  additional  unknowns  u.u.,~t  ' (3u  /3x  ) 

1  J  ij  i  j 

and  u^9.  Models  proposed  so  far  to  evaluate  these  unknowns  have  them 
coupled  to  the  mean  quantities  through  either  algebraic  or  differential 
equations.  Some  are  based  on  empirical  relation  and  others  on 
postulations . 


In  1877,  Boussinesq  [2]  proposed  the  concept  of  eddy  viscosity  which 
assumes  that,  in  analogy  to  the  viscous  stresses  in  laminar  flows, 
turbulent  stresses  are  proportional  to  the  mean  velocity  gradients .  For 
general  flow  situations,  it  is  expressed  as 


_  3U.  au. 

-u.u.  =  v  r — -  +  — li  -  -vs 
i  j  tl3x.  3x .  3 KOij 

J  i  J 


(1.7) 
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vt  is  the  turbulent  or  eddy  viscosity  which,  unlike  molecular  viscosity, 

is  not  a  fluid  property  but  depends  on  the  state  of  turbulence,  k 

represents  the  kinetic  energy  of  the  fluctuating  motion  or  u  u  7 T.  The 
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above  expression,  however,  does  not  close  the  problem  of  turbulent  flows 
as  and  k  are  still  unknowns.  In  1925,  Prandtl  [3]  proposed  a 
turbulence  model  called  the  'mixing  length'  model.  This  model  provides  a 
relation  between  the  eddy  viscosity,  a  length  scale  L  characterizing  the 
size  of  turbulent  eddies  and  a  suitable  velocity  scale,  V.  Thus 

vt  «  V*L 

Both  the  turbulent  velocity  scale,  V,  and  the  mixing  length  scale,  L, 
could  be  reasonably  approximated  for  many  flows.  However,  for  such 
flows,  empirical  constants  were  needed  to  prescribe  a  length  scale.  In 
most  of  these  flows,  the  constants  were  obtained  by  fitting  the 
calculated  results  to  experimental  data  of  a  particular  flow  under 
study.  These  mixing  length  model  constants  were  found  [4]  to  vary  often 
from  one  flow  to  another.  Consequently,  the  mixing  length  turbulence 
model  is  successful  only  in  predicting  turbulent  flows  in  similar 
geometry  and  flow  conditions  but  lacks  the  universality  and 
predictability  when  the  turbulent  flow  and  geometry  conditions  are 
different.  Other  models  [5,6],  similar  to  the  mixing  length  model,  were 
shown  to  have  success  in  a  given  flow  but  lacked  generality  when  flow 
conditions  and  configuration  changed. 
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To  overcome  the  lack  of  predictability  and  generality,  several  more 
complex  models  [7,8]  were  developed  during  the  1940's  and  1950's  which 
employed  differential  transport  equations  for  the  turbulent  quantities. 
However,  these  equations  could  not  be  solved  directly  as  there  were 
mathematical  difficulties  involved  and  numerical  techniques  and  fast 
computers  were  not  available.  Alternatively,  the  governing  partial 
differential  equations  for  turbulent  flows  were  often  integrated  and 
reduced  to  ordinary  differential  equations.  These  integral  methods 
assumed  some  shape  of  mean  profile  and  used  some  empirical  relations  for 
global  behavior  of  turbulence.  They  lacked  flexibility  since  the  assumed 
profile  must  be  approximately  the  same  in  the  flow  field  and  could  not 
be  applied  for  different  flows. 


Advances  in  computational  facilities  and  numerical  methods  during  the 
late  1960  s  and  1970's  led  to  the  use  of  more  advanced  models  which 
solve  complete  partial  differential  equations  for  both  mean  flow  and 
turbulent  quantities.  One  of  these  models  which  solves  the  differential 
equation  for  k,  the  kinetic  energy,  is  called  the  one-equation  model  as 
opposed  to  the  zero-equation  model  where  no  differential  equations  are 
solved  for  turbulent  quantities.  With  the  kinetic  energy  known,  the 
eddy  viscosity  can  be  written  as 


V  =  C  k^L 
t  y 


where  represents  a  velocity  scale,  L  the 
empirical  constant.  The  equation  for  k  is 


length  scale  and 


(1.8) 


an 


10 


Dk  _  3 _  \  3k_ 

Dt  3x .  La,  3x . 

J  k  j 


u  ,u . 
1  J 


au. 

_ i 

3x . 
J 


1  *  5 


Cl. 9) 


which  is  derived  from  the  governing  equation  of  fluctuating  turbulent 

motion.  Details  of  the  derivation  are  given  later.  Here.  and  a  re 

D  k 

empirical  constants.  This  one-equation  model  is  not  complete  unless  the 
length  scale  L  is  specified.  In  most  cases,  L  is  a  variable  and  is 
obtained  from  simple  empirical  relations  similar  to  those  for  the  mixing 
layer . 


Since  one-equation  models  [9,10]  account  for  the  convective  and 
diffusive  transport  of  the  turbulent  kinetic  energy,  they  are  superior 
to  the  mixing  length  models  in  flows  where  the  transport  mechanism  is 
important.  Some  examples  are  non-equilibrium  boundary  layers  with 
rapidly  changing  free-stream  conditions,  boundary  layers  with  free- 
stream  turbulence  and  recirculating  flows.  However,  in  many  flows  it  is 
difficult  to  specify  the  length  scale  empirically.  The  logical  extension 
of  the  turbulence  modelling  is  that  the  length  scale  be  obtained  from  a 
differential  transport  equation. 


Models  which  solve  differential  equations  for  both  turbulent  velocity 
scale  or  turbulent  kinetic  energy  and  length  scale  are  known  as  two- 
equation  models.  Several  different  models  [4,11]  have  been  proposed 
which,  in  addition  to  the  equation  for  k,  solve  an  equation  of  the  form 
k^11  instead  of  L.  The  most  popular  one  is  the  one  suggested  by  Jones 
and  Launder  [11]  which  has  m=1.5  and  n=-l.  The  term  k1‘5L*1  which 
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appears  in  the  last  term  of  equation  (1.9),  has  a  physical  significance 
as  it  has  the  same  dimension  as  e,  the  dissipation  of  turbulent  energy. 
The  dissipation  function  of  turbulent  kinetic  energy,  z  or 
v(3u^/3x.) (3u./3x.)  can  be  derived  and  modelled  as 

^  J  J 


Dz 


k2  3s 


g  _  3U.  a 

Dt  3x  J-Cs  s  3xJ  Csl  k  UiUj  3x7  "  Cz2  k 
J  J  J 


(1.10) 


Details  of  the  derivation  of  equation  (1.10)  are  given  later.  Here,  C 

anc^  ^£2  are  emParaca^  constants .  The  k-s  model  with  eddy  viscosity 

from  equation  (1.7)  now  requires  six  emipirical  constants  C  .  a  C 

y  k’  D’ 

V  C£1  “d  °c2' 


This  k-E  model  has  been  used  in  the  calculation  of  boundary  layer 
type  of  flows  as  well  as  recirculating  flows.  The  model  now  can  predict 
large  number  of  different  flow  configurations  and  conditions  and  is 
certainly  more  general  than  the  mixing  length  turbulence  model.  Though 
this  model  has  a  wider  range  of  application  in  the  past  fifteen  years, 
it  still  lacks  universality  as  the  coefficients  need  to  be  adjusted  from 
one  flow  to  another.  As  an  example,  the  constant  Cg2  in  the  e -equation 
has  a  value  between  1.90  and  1.92.  Using  this  value  of  C£?,  a  reasonably 
good  prediction  of  plane  jet  flow  can  be  made.  However,  if  the  value  of 
this  constant  is  slightly  outside  this  range,  the  solution  becomes 
sensitive  to  the  constant  and  does  not  converge.  Furthermore,  the  value 
^ £2  between  1.90  and  1.92  which  gives  good  prediction  of  plane  jet 
flow  cannot  be  used  for  a  round  jet  since  it  produces  a  30%  error  in  the 
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spread  of  turbulent  round  jet.  For  a  round  jet,  the  value  of  is 

found  not  to  be  a  constant  and  is  changed  [4]  to  1 . 9 2* (1-0 . 035G)  where 


3U 


(1.11) 


Another  problem  is  that,  if  these  modified  k-E  model  equations  (1-7), 
(1.8),  (1.9),  (1.10)  and  (1.11)  were  used  for  the  calculation  of  plane 
wake  flow,  there  is  a  30%  under-prediction  in  the  growth  or  spread  of 
the  wake.  This  difficulty  is  further  taken  care  of  by  making  the 
constant  in  equation  (1.7)  a  function  of  P/e  [4]  where  P  is  the 


production  of  turbulent  kinetic  energy  -u^u^ (3U^/3x^ )  and  z  is  the 


dissipation  of  this  energy,  v(3u^/3x^  )  (3u^/3x^. ) 


It  should  be  remarked  here  that  these  difficulties  are  mainly  dealing 
with  the  generality  or  universality  of  the  model.  In  general,  the  k-E 
model  has  achieved  a  level  of  predictability  which  mixing  length  or  one- 
equation  turbulence  models  could  not.  In  order  to  advance  the 
predictability  of  turbulent  flow  motion  further  improvement  in 
turbulence  modelling  must  be  made.  This  motivates  the  present 
investigation. 


1.3  Scope  of  the  Present  Work 

In  this  investigation,  a  fundamental  change  in  turbulence  modelling 
is  made,  that  is,  to  introduce  the  two  scale  concept,  one  based  on  (k,£) 
scale  and  the  other  (e,v)  scale.  In  the  present  investigation,  k  and  z 
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are  used  to  scale  the  turbulent  phenomenon  dominated  by  large  scale 
motion  such  as  diffusion  term  while  the  physical  process  associated  with 
the  dissipation  of  turbulent  kinetic  energy  is  modelled  using  e  and  v  as 
the  basic  parameters,  which  is  known  as  Kolmogorov  scale.  The  Kolmogorov 
scale  which  is  known  since  1925  is  more  closely  related  to  small  eddy 
motion  and  has  not  been  incorporated  in  the  turbulence  modelling  so  far. 
However,  in  the  present  investigation,  this  scale  is  used.  The  new 
turbulence  model  based  on  both  (k,e)  and  (e,v)  scale  is  called  the  two- 
scale  turbulence  model. 

In  Chapter  II,  a  description  of  the  physics  of  turbulence  and  the 
theory  behind  the  use  of  the  two-scale  model  is  given.  Then,  the 
detailed  derivation  of  the  two-scale  turbulence  model  is  shown.  Chapter 
III  gives  the  governing  equations  for  buoyant  flows.  Chapter  IV  contains 
a  review  and  collection  of  experimental  data  for  free  shear  flows.  In 
Chapters  V  and  VI  the  prediction  of  several  free  shear  flows  is  shown. 
Chapter  VII  shows  the  calculations  for  separated  flows.  Finally,  chapter 
VIII  contains  several  important  observations  about  the  model  and 
possible  areas  of  further  work  regarding  multiple  scale  modelling. 
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CHAPTER  II 

TWO -SCALE  SECOND  ORDER  TURBULENCE 
MODEL  FOR  INCOMPRESSIBLE  FLOWS 

This  chapter  gives  a  detailed  derivation  of  the  two-scale  k-e  model 
for  incompressible  turbulent  flows.  The  complete  set  of  governing 
equations  are  presented  which  are  then  modelled  based  on  a  set  of 
turbulent  postulations. 


2.1  Governing  equations 

The  governing  equations  for  incompressible  turbulent  flow  are  the 
averaged  Navier-Stokes  equations,  namely,  the  continuity  equation,  the 
momentum  equation  and  the  energy  equation.  They  are  also  known  as  the 
Reynolds  equations  since  it  was  Reynolds  [1]  who  first  used  the 
averaging  technique.  For  a  short  time  or  ensemble  average,  the  average 

ic 

value  of  an  instantaneous  quantity  0  at  a  time  t  can  be  defined  as 


Z  **(*»“) 
n=l 

th. 

where  n  denotes  the  n  measurement  of  a  total  of  N  experiments .  In 
cartesian  tensor  notations,  the  continuity  equation  is 


(2.1) 
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The  momentum  equation  is 


DU. 
P _ i 

Dt 


+ 


(2.2) 


where  x..  and  x..  are  the  laminar  and  turbulent  stresses,  G.  is  the 
1J  i 

body  force  and  P  is  the  pressure.  The  stresses  x  _  and  x  t  are  given  by 
the  relations 


x 


au. 


and 


pu  .u . 
i  J 


The  term  -pu^u^,  known  as  Reynolds  stress,  is  a  result  of  averaging  the 
convective  acceleration.  It  is  generally  regarded  as  a  turbulent  stress 
in  analogy  with  viscous  stress,  and  is  unknown.  The  energy  equation, 
which  too  has  additional  unknown  quantities,  is  given  by 


DT 

pCDt 


au. 

T  .  .  1 

1JaT 
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9x . 
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+  * 


(2.3) 


where  the  laminar  heat  flux  q_^  and  the  turbulent  heat  flux  q.^  are  given 
by  the  relations 


q.  -  -k—  and  q.  =  -pcu.0 

l  8x .  r  i 

l 

0  is  the  viscous  dissipation  due  to  the  velocity  fluctuation  and  is 
expressed  as 
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3u  .  3u  .  3u  . 

=  u(aT  + 

J  1  J 


In  the  above  five  equations  there  are  fifteen  unknowns,  namely,  ,  P, 
T,  uiuj  >  u^9  and  0.  Hence,  it  is  necessary  to  obtain  equations  fcr 
u^Uj  ,  uTH  and  0  to  complete  the  turbulence  closure  problem. 

Equations  for  fluctuating  velocity,  u^,  and  fluctuating 
temperature, 8 ,  are  obtained  by  subtracting  the  above  averaged  equations 
from  the  original  Navier-Stokes  equations.  This  gives  the  momenrura 
equation  denoted  by  (ra^)  for  the  fluctuating  velocity  component,  u^, 


Du .  3u .  3u . 

l  l 


8uiut  .  1  .  3‘Ui 

p  ax,  “aiyf 


(m,)(2.4) 


and  the  energy  equation  denoted  by  (0)  for  fluctuating  temperature,  8, 


D8  .  3T 
Dt  +  nl  W. 


+  U 
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pc 
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where 


3u .  3u .  3u .  3U .  3u .  3u .  3u . 

''  =  Tii  ST  +  +  a^'tT  +  «<*r  + 

J  J  1  J  J  X  J 


From  equation  (2.4),  the  equation  for  u^u^  is  obtained  using  the 

relation  f(m . )u .+(m . )u . ] .  This  results  in 
1  J  J  1 
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(2.6) 


In  the  above  equation,  the  first  term  on  the  right  hand  side  represents 

both  the  molecular  and  turbulent  diffusion  of  the  stress  u  u  The  next 

i  J 

term  is  the  product  of  the  Reynolds  stress  and  the  strain  rate  which 
represents  the  interaction  between  fluctuating  component  and  mean  flow. 
It  is  often  called  the  production.  The  third  terra  is  the  dissipation. 
The  last  term  in  this  equation  represents  the  correlation  between 
pressure  and  fluctuating  velocity  gradients.  It  is  also  called  the 
pressure-strain  term  or  the  redistribution  term.  The  above  equation  can 
be  contracted  to  get  the  equation  for  turbulent  kinetic  energy  k  or 
u£u^/2  by  summing  i=j  and  dividing  it  by  2.  This  gives  with  s  = 
vl(3u./3xjl)Caui/3xjl)] 
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(2.7) 


where  the  term  on  the  left  side  represents  the  time  rate  change  of 

turbulent  kinetic  energy  following  the  mean  convection  U  .  The  first 

i 

term  on  the  right  side  is  the  diffusion  of  k.  The  second  and  third  terms 
are  the  production  and  dissipation  of  the  turbulent  kinetic  energy.  The 
dissipation  term,  z ,  represents  the  rate  of  dissipation  of  turbulent 
kinetic  energy  and  is  an  unknown  in  the  above  equation.  It  should  be 
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remarked  that  the  dissipation  term  z  appears  naturally  in  the  k- 
equation.  The  variation  of  z  in  the  flow  field  has  an  important  bearing 
of  the  distribution  of  the  turbulent  kinetic  energy.  Thus  z  is  an 
important  turbulent  transport  property.  The  differential  equation  for  z 
is  derived  from  the  (m^)  equation  by  using  the  relation 
2v[  3(nu  )/3x^]  [3u^/3x^]  .  This  gives 
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(2.8) 


It  should  be  noted  that  although  the  above  equation  is  exact  but  every 
term  on  the  right  side  other  than  the  viscous  diffusion  v(3e/3x^),  is  an 
additional  unknown  quantity.  The  first  term  on  the  right  side  is  the 
diffusion  of  e  while  the  second  and  third  terms  represent  the  production 
of  e.  The  last  two  terms  are  often  called  the  destruction  of  z.  The 
modelling  of  these  terms  will  be  done  in  the  next  section. 


Finally,  the  u^0 -equation  is  obtained  from  equations  (2.4)  and  (2.5) 
by  using  the  relation  [QCmp+u  (8)]  which  results  in 
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where 
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=  T  .  . 


3u .  Bu .  3u .  3U .  Bu .  Bu  .  Bu  . 

ij  3x  +  Ml3x  +  3x.  hx.  +  P^3x.  +  SxJsxT 
J  J  1  J  j  i  j 


In  this  equation,  the  terms  on  the  right  side  are  diffusion  of  u , 8  the 

i 

production  of  u~TT^  the  dissipation,  the  pressure-temperature  correlation 

and  the  frictional  heating  terms  respectively.  The  unknown  0 1 u .  in  the 

uTS - equat ion  represents  the  frictional  -heating  generated  by  the 

fluctuating  component  and  is  usually  considered  to  be  smaller  than  the 

frictional  heating  generated  by  the  mean  flow  motion  t  (3U  /3x  ) 

ij  i  j 

Hence,  it  is  often  omitted  in  the  mean  energy  equation.  It  should  also 

be  noted  here  that  a  part  of  the  mean  energy  equation  (1.6)  yfBu  /3x  + 

_ 1  j 

3u^/3x^] I3u^/3x^]  is  equal  to  z  which  is  derived  in  equation  (2.8). 


The  four  transport  equations  (2.6)  to  (2.8)  derived  above  have 
several  unknown  terms  on  the  right  side  most  of  which  need  to  be 
modelled.  This  is  discussed  in  the  following  section. 


2.2  Concept  of  Two  Turbulent  Scales 

Before  attempting  to  model  these  equations,  a  brief  discussion  of 
turbulent  flow  structure  is  done  and  the  concept  of  the  two  turbulent 
scales  is  introduced.  In  order  to  visualize  the  existence  of  two 
significantly  different  turbulent  scales  in  a  turbulent  flow,  it  is 
instructive  to  consider  a  turbulent  correlation  function  R  (x;r)  for 

ij 

velocity  fluctuation,  which  is  defined  as 


Ri  (x;r)  =  u  (x)  u.(x+r) 
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where  u^(x)  is  the  instantaneous  value  of  the  i^  component  of  the 
fluctuating  velocity  at  the  point  of  the  position  vector  x  and  u^.(x+r) 
the  component  of  the  fluctuating  velocity  at  (x+r) .  The  average, 

with  a  bar  over  u^uj  maY  be  considered  either  a  time  average  or  an 
ensemble  average.  If  r=0  and  i= j ,  the  one  point  correlation  R^CxjO)  is 
the  Reynolds  normal  stress  in  the  it^1  direction.  The  correlation 
R^(x,0)  includes  all  possible  turbulent  eddy  sizes  at  the  position  x. 

It  is  difficult  to  differentiate  the  scale  that  is  significant  in 
carrying  out  a  turbulent  process.  One  way  to  examine  the  behaviour  of 
each  turbulent  eddy  is  to  consider  a  spectral  analysis  of  the 


correlation  R..(x,0),  i.e. 
ij 


Vk) 


■  cfc’f 


R^(x;r)  exp(-ik.r)  dr 


where  (k.r)  is  the  wave  number  vector,  k  dot  the  position  vector  at  r 
distance  from  x.  The  wave  number  vector  may  be  written  as 


k  =  ki  +  kj+kk 
x  yJ  z 

The  component  wavenumber,  k^,  is  related  to  the  fluctuating  frequency  n^ 
and  the  wavelength  1,  of  an  eddy  in  the  x^  direction  by 


o  2irn . 

k  _  2it  _ i 

l  X.  “  U. 

1  x 


In  fact  0 . . (k)  is  the  Fourier  transformation  of  R..(x;r).  The  inverse 
ij  ij 

Fourier  transformation  for  recovering  R„(xfr)  thus  becomes 
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R.  .(x;r)  -  f 0 . . (k)  exp(ik.r)dk 
1J  „  1J 


The  reason  for  examining  the  spectral  distribution  is  that  the 
transform  is  simply  a  method  of  representing  the  complex  random  wave 
form  of  turbulent  eddy  motion  associated  with  R„  by  what  is  equivalent 
to  a  sum  of  sine  or  cosine  waves  of  various  amplitude  or  frequencies. 

The  total  sum  of  all  sine  and  cosine  waves  is  equivalent  to  the  original 
wave  form  of  R_(x;'r).  Thus,  one  may  think  of  0  (£)  as  a  fluctuating 

intensity  of  R^(x;r)  at  a  wave  number,  k^,  or  frequency  n^.  If  the 
fluctuating  intensity  is  large  at  a  particular  range  of  wave  numbers,  it 
means  that  the  physical  process  of  the  turbulent  phenomenon  is 
intimately  related  to  this  range  of  wave  number. 

For  the  present  analysis,  the  energy  spectrum  of  a  steady  isotropic 
flow  behind  a  wind  tunnel  grid  at  r=0  is  considered.  Then 


QC 

R  (x;0)  (k)dt 

—  At  ^ 


The  energy  spectral  is  a  function  of  the  wave  vector  k  or  of  a 

given  point  at  k  in  wave  space.  An  integrated  energy  spectrum  E^Ck) 
which  is  a  function  of  a  scalar  variable  k  can  be  obtained  by 
integrating  the  energy  spectrum  0_(k)  over  a  spherical  surface  of 
radius  k=|k|  or 

Eij(k)  =  j*ijOOds(k) 
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Here,  ds(k)  is  an  element  on  the  surface  of  the  sphere  of  radius  k. 
E„(k)  thus  may  be  taken  as  the  energy  contribution  from  the  eddy  size 
with  wave  number  k  to  the  u^u^  correlation.  The  energy  spectrum  function 
of  turbulent  kinetic  energy  in  the  wave  space  is 


E(k)  =  fE..(k) 


The  total  kinetic  energy  of  the  turbulent  flow  is  then 


u  .  u . 
1  1 

2 


*0 

E  (k)dk 
0 


In  particular,  for  isotropic  flow  the  relation  is 

oO 

f E(k)dk  =  |u2 

The  spectrum  equation  of  turbulent  kinetic  energy  equation  for  isotropic 
turbulence  can  be  written  [12]  as 


3E(k) 

3t 


=  T(k)  -  D(k) 


where  T(k)  is  associated  with  the  transfer  of  energy  between  wave 
numbers  or  eddy  sizes.  Its  integral  over  all  wave  numbers  is  zero.  It 
can  thus  be  defined  by  a  different  transfer  function 

(k 

S(k)  =  -  I  T(k)dk 

-'O 
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which  is  the  total  energy  transferee!  from  eddies  in  the  range  from  0  to 
k  to  those  in  the  range  greater  than  k.  In  other  words,  S(k)  is  the  flux 
turbulent  kinetic  energy  from  a  spherical  volume  of  radius  equal  to 
wave  number  k.  D(k)  is  the  rate  of  dissipation  of  turbulent  kinetic 
energy  at  the  wave  number  k  and  is  equal  to 

D(k)  =  2vk2E(k) 


Figure  2.1  shows  the  schematic  energy  spectrum  E(k,t)  and  the 
dissipation  spectrum  D(k,t)  for  an  isoptropic  flow.  The  solid  line  shows 
a  typical  energy  spectrum  and  the  dashed  line  the  dissipation  of 
turbulent  kinetic  energy.  Figure  2.2  gives  the  measured  energy  spectrum 
and  the  dissipation  [  12-14  ]  in  log-log  scale  for  a  steady  flow  behind 
a  square  grid  screen  with  spacing  of  M  in  a  wind  tunnel.  Here,  the 
dimensionless  wavenumber  k  is  defined  as  2Trnti/U  with  n  the  frequency  of 
a  fluctuating  component  in  turbulent  flow,  U  the  mean  flow  velocity  and 
Tl  is  the  Kolmogorov  length  scale  or  (v?/e)^.  t  is  a  dimensionless  time 
or  the  real  time  normalized  by  a  characteristic  time  M/U.  In  figure  2.2, 
the  Reynolds  number  Rex  is  UX  /v  where  X  is  Taylor's  microscale  [12]. 
The  wavenumber,  k,  may  be  considered  to  be  inversely  proportional  to  the 
size  of  the  eddies.  In  other  words,  the  larger  the  size  of  the  eddy,  the 
smaller  is  its  wavenumber.  From  figure  2.2,  it  can  be  seen  that  the 
measured  energy  and  dissipation  spectra  are  quite  different  and  can  be 
associated  with  different  wavenumbers.  For  instance,  a  wavenumber 
characterized  by  k^,  in  the  order  of  10  ^  at  Reynolds  number  Re^  of  540 
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may  be  considered  to  be  associated  with  the  size  of  the  small  eddies 
that  provide  the  main  contribution  to  the  dissipation  of  turbulent 
kinetic  energy.  This  value  k^  roughly  corresponds  to  the  maximum  value 
of  the  dissipation  curve.  Similarly,  there  is  a  range  of  spectrum  which 

corresponds  to  the  energy  containing  large  eddies.  A  wavenumber 

-4 

characterized  by  k  ,  in  the  order  of  10  at  Re.  of  540  may  be 

e  a 

considered  to  associate  with  this  range  which  corresponds  to  the  peak  of 
the  energy  curve. 


It  has  been  shown  both  experimentally  by  Frieche  et  al.  [13]  and 
theoretically  by  Driscoll  and  Kennedy  [14]  that  these  energy  and 
dissipation  spectra  change  with  Reynolds  number.  As  given  in  figure 
2.2,  an  increase  in  the  Reynolds  number  causes  the  peaks  of  the  energy 
and  dissipation  curves  to  separate  further  away. 


In  most  of  the  spectral  analysis,  a  turbulent  Reynolds  number  is 

associated  with  the  wavenumber,  k  .  It  has  been  shown  [12]  that 

e  L  J 


k~  =  l  =  T^Re.X 
k  e  15  X  g 
e 


where  A  is  a  constant  and  i is  the  length  of  the  eddy  corresponding  to 

the  wavenumber,  k^.  Re^  is  the  Reynolds  number  based  on  Taylor 

microscale,  X  or  UX  /v.  The  Taylor  microscale  is  a  length  scale 
S  S 

associated  with  the  curvature  of  the  spatial  velocity  autocorrelations 
[15]  and  is  related  to  the  dissipation  t  by  the  expression  [12] 


i  -  u 

£  =  lDVX 


»  2 
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where  u'  is  a  velocity  fluctuation. 

Driscoll  and  Kennedy  [14]  obtained  the  energy  and  dissipation  spectra 

for  Re^  ranging  from  13  to  540  as  shown  in  figure  2.2  The  dimensionless 

wavenumber,  k,  is  defined  as  2irnr|/U  where  ^  is  the  Kolmogorov  length 

scale  or  (v3/e)^.  The  energy  spectra  shows  that  when  Re.  increases  k  n 

a  e 

decreases.  For  a  value  of  Re. =13,  the  peak  wave  number  k  is  about  0.01 

*  e 

whereas  for  Re^=540,  it  is  0,0001.  Hence,  it  can  be  said  that  the 
structure  of  turbulence  is  dependent  on  Reynolds  number,  whether  it  is 
the  turbulent  Reynolds  number  or  the  mean  Reynolds  number. 

The  quantity  E(k)  [12] ,  used  in  Figures  2.1  and  2.2,  is  defined  as 

E(k)  =  27Tk2E.  .(k) 

li 

where  E^Ck)  is  the  Fourier  Transformation  of  the  correlation  tensor 

u  .u .  or 
l  l 

EiiOO  =  (' uiuiexp(-ik."r)  dr 

Thus,  the  total  energy  contained  by  all  the  eddies  is  l^ll^.i.e. 

E  (k)dk  = 

Therefore,  figure  2.1  shows,  conceptually,  two  distinguishing  features 
of  turbulence  when  one  examines  the  turbulent  spectra  or  turbulent 
eddies.  The  solid  line  gives  the  energy  spectra  from  which  it  can  be 
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Figure  2.1.  Energy  and  dissipation 
spectrum  of  an  isotropic  flow 
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Figure  2.2.  Energy  and  dissipation 
spectra  for  various  Re. 
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seen  that  the  range  of  eddies  containing  most  of  the  energy  are  large  in 
size  (or  lower  in  wavenumber  range)  and  is  comparable  to  the  width  of 
the  flow.  They  transfer  their  energy  to  smaller  eddies.  It  is  in  this 
range  of  smaller  eddies  where  most  of  the  dissipation  of  turbulent 
energy  occurs.  The  larger  the  Reynolds  number,  the  smaller  is  the  eddy 
size.  These  properties  of  turbulent  flows  are  obtained  by  experimental 
measurements  and  not  by  any  postulation.  Hence,  it  seems  natural  to 
consider  different  scales  for  the  modelling  of  the  k  and  z  equations 
(2.7)  and  (2.8).  The  measurements  of  Frieche  et  al.  [13]  reveal  that 
large  eddies  possess  most  of  the  turbulent  kinetic  energy  in  the  flow 
and  do  not  play  any  significant  role  in  the  dissipation  of  turbulent 
kinetic  energy.  On  the  other  hand,  Kolmogorov  [12]  found  that  small  eddy 
characteristics  are  functions  of  (s,v).  In  the  medium  range  of  eddy 
size,  a  process  described  as  the  transfer  function  T(k,t)  derived  from 
convection  terms  of  the  k-equation  (2.7)  provides  a  mechanism  to 
transfer  the  turbulent  kinetic  energy  possessed  by  large  eddies  to  small 
eddies  before  it  is  consumed  by  the  viscous  dissipation  and  turned  into 
thermal  energy.  This  distinct  difference  in  the  behavior  of  turbulence 
different  wave  number  was  known  for  sometime.  However,  it  has  not  vet 
been  incorporated  in  most  of  the  turbulence  models .  The  existing  models 
characterize  the  velocity,  length  and  time  scales  for  turbulent  flows 
based  on  k  and  e.  However,  in  any  turbulent  flow,  it  is  the  larger 
eddies  which  cascade  to  become  smaller  eddies  through  inertial 
interaction,  thereby  transferring  energy  to  the  smaller  eddies.  At  the 
same  time,  viscosity  effects  and,  with  them,  dissipation  become  more  and 
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more  important  for  the  smaller  eddies  as  shown  in  figure  2.1.  For  a 
certain  range  of  these  small  eddies ,  it  can  he  shown  that  turbulence  is 
in  statistical  equilibrium.  This  is  the  range  in  which  viscosity  can  be 
effective  in  smoothing  out  velocity  fluctuations.  The  generation  of 
these  small  scale  fluctuations  is  made  possible  due  to  the  nonlinear 
terms  in  the  equations  of  motion.  On  the  other  hand  the  viscous  action 
prevents  the  generation  of  infinitely  small  scales  of  fluctuating  motion 
by  dissipating  turbulent  kinetic  energy  into  heat.  One  may  consider  that 
at  large  Reynolds  numbers,  the  relative  magnitude  of  viscous  force 
compared  to  inertia  force  is  so  small  that  viscous  effects  in  a  flow 
tend  to  become  vanishingly  small.  However,  Townsend  [15]  reasoned  that 
the  nonlinear  terms  in  the  Navier-Stokes  equations  counteract  this 
effect  by  generating  motion  at  scales  small  enough  to  be  affected  by 
viscosity.  In  other  words,  as  soon  as  the  scale  of  the  flow  field 
becomes  so  large  that  viscosity  effects  could  be  neglected,  the  flow 
creates  small  scale  motion  thereby  keeping  viscosity  effects  and,  in 
particular,  dissipation  rates  at  a  finite  level. 

At  these  small  scales,  turbulent  motions  are  statistically 
independent  of  the  relatively  slow  large  scale  turbulence  and  of  the 
mean  flow.  Hence,  as  Kolmogorov  reasoned,  the  character  of  turbulence  in 
this  range  is  determined  by  s,  the  rate  of  dissipation  of  k  and  the 
viscosity  v.  These  considerations  led  Kolmogorov  to  make  the  following 
hypothesis : 
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At  suf ff icient ly  high  Reynolds  numbers  there  is  a  range  of  high  wave 
number  where  the  turbulence  is  statistically  in  equilibrium  and  uniquely 
determined  by  the  parameters  z  and  v.  This  state  of  equilibrium  is 
universal 1 . 

Using  these  two  parameters,  z  and  v,  velocity,  length  and  time  scales 
for  small  eddy  motion  can  be  characterized  by 

v  =  (ve)^;  T)  =  (v3/04  and  t  =  T|/v  =  (v/s)^ 

which  can  be  obtained  by  dimensional  analysis  of  v  and  z .  On  the  other 
hand,  in  the  large  turbulent  eddies,  the  turbulent  kinetic  energy,  k,  is 
important  since  these  large  eddies  are  responsible  for  carrying 
turbulent  energy  and  extract  energy  from  the  flow  motion  to  sustain 
turbulence.  Therefore,  the  character  of  turbulence  in  the  large  eddy 
range  is  determined  by  e,  the  rate  of  dissipation  of  k,  and  the 
turbulent  kinetic  energy,  k,  itself.  Using  these  two  parameters,  z  and 
k,  the  velocity,  length  and  time  scales  for  large  eddy  motion  can  be 
characterized  as 

u  =  k^;  l  =  kl's /z  and  t  =  i/u  =  k/e 

Though  the  above  analysis  was  done  for  isotropic  flow,  which  is  not 
the  case  in  many  practical  situations,  it  has  been  shown  [12] 
experimentally  that  the  fine  structure  of  nonisotropic  turbulent  flows 
is  almost  isotropic  (local  isotropy).  This  is,  however,  not  true  for  all 
experimental  results.  Nevertheless,  many  qualitative  features  of 
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isotropic  turbulence,  particularly  the  distribution  of  two  turbulence 
scales,  apply  to  phenomena  in  actual  turbulence.  Measurements  of 
Kolmogorov  fine-scale  turbulence  structure  in  various  flows  shows  that 
differences  between  results  are  often  sufficiently  small  to  be 
negligible  in  the  first  approximation. 

Several  investigators  [16,17]  have  mentioned  in  the  past  that  it  is 
the  t -equation  [equation  (2.8)]  which  needs  to  be  carefully  studied. 

This  is  because  of  the  complexity  and  difficulty  in  modelling  the  £- 
equation.  The  physical  meaning  of  the  different  correlations  among  all 
sizes  of  eddies  and  fluctuating  quantities  is  sometimes  difficult  to 
understand.  As  an  example,  the  production  term  containing  the  second 
derivative  of  the  mean  velocity  U  in  equation  (2.8)  for  e  is  neglected 
invariably  by  most  investigators.  The  reason  for  this  is  that  this  term 
is  assumed  to  be  much  smaller  than  some  of  the  other  terms  in  this 
equation.  However,  the  physical  significance  of  this  term  is  still  not 
clear.  Therefore,  due  to  lack  of  information  about  such  terms  the  e- 
equation  needs  to  be  further  investigated  in  order  to  improve  the 
accuracy  and  prediction  capability  of  the  model  as  well  as  making  it 
more  general. 

The  concept  of  using  different  time  scales  was  first  proposed  by 
Lumley  [17]  in  1975.  He  suggested  that  each  term  in  both  the  k  and  e 
equations  be  modelled  either  by  using  the  (k,s)  scales  or  the  (k,s,v) 
scales.  However,  in  the  final  form  of  the  modelled  £ -equation  suggested 
by  Lumley,  the  scale  containing  v  was  neglected.  Another  approach 


32 


considering  multiple  scales  was  made  by  Hanjalic  et.  al  [16].  They  used 
two  different  time  scales  by  dividing  the  whole  energy  spectrum  into  two 
parts  --  the  energy  containing  eddies  and  the  dissipating  eddies.  For 
each  region,  a  separate  time  scale  is  used  to  model  the  k  and  z 
equations.  Results  were  obtained  for  several  thin  shear  flows  which  show 
an  improvement  in  the  level  of  agreement  wTith  experiments  over  that 
obtained  with  models  employing  only  one  time  scale.  The  authors 
suggested  that  by  dividing  the  spectrum  into  more  number  of  parts  and 
solving  the  two  equations  in  each  region,  a  further  improvement  in  the 
result  could  be  obtained  though  the  computational  time  would 
considerably  increase.  However,  the  authors  did  not  present  the 
results . 

In  the  present  investigation  of  turbulence  modelling,  the  two-scale 
concept  is  employed.  The  two  scales  are  the  large  eddy  or  energy 
containing  scale  based  on  k  and  z  and  the  small  eddy  or  energy 
dissipating  scale  based  on  v  and  z.  The  two-scale  concept  is  applied  to 
all  turbulent  transport  equations  whenever  it  applies. 

2.3  Turbulence  modelling 

Before  modeLling  the  transport  equations,  the  postulations  of 
turbulent  flow  are  listed  below.  These  postulations  are  made  by  various 
models  and  summarized  by  Chen  [18]. 

1.  Navier-Stokes  equations  are  valid  in  describing  turbulent  motion. 
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Turbulent  diffusion  of  a  turbulent  transport  quantity  (u.u.,  k, 


3. 

4. 


£,  uTo)  is  proportional  to  its  gradient. 
Small  eddies  are  isotropic. 


All  turbulent  quantities  are  functions  of  u.u.,  k,  e .  u,6.  U  ,  P, 

i  j  i  i*  * 

T ,  p ,  v  and  a . 

5.  The  model  equations  should  be  consistent  with  respect  to 
symmetry,  invariance,  permutation  and  physical  conservation  laws 
imposed  on  the  original  equations. 

6.  Turbulent  scales  are  functions  of  k,  e  and  v.  Large  eddy  scales 

based  on  (k,e)  are  [u]  =  k2 ,  [1]  =  k1'5/^,  [t]  =  k/e  and  small 

eddy  scale  based  on  (v,e)  are  [u]  =  (ve)^’  ^  ~  ^ 

(v/e)^. 

7 •  Turbulent  constants  in  the  model  are  determined  from  experiments . 

The  two-scale  turbulent  flow  model  is  now  derived  in  the  following 
section.  Both  (k,s)  and  (v,e)  scales  are  used  in  the  modelling  of  the  e- 
equation.  As  for  the  modelling  of  the  U7u\  and  k  equations,  the  large 
eddy  scale  (k,e)  is  used  for  the  reason  that  the  large  eddies  which 
contain  most  of  the  turbulent  kinetic  energy  are  also  responsible  for 
turbulent  diffusion  and  pressure-strain  interaction.  Further  details 
are  presented  below. 


2.3.1  Modelling  of  UTu*.  and  k  equations 
The  turbulent  diffusion  term  of  equation  (2.6)  is  modelled  based  on 
postulate  2  that  the  diffusion  of~u“u~  is  proportional  to  its  gradient 
or 
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,2  9u.u. 


-  u.u.u  +  -(6..U.  +  6..u.)  =  C.  [  —  ]  — 
1  j  l  p  j£  i  i£  j  kLt  J  3x 


In  order  to  keep  the  dimensions  consistent,  a  quantity  with  a  scale  of 
[l2/t]  is  needed  to  complete  the  model.  From  dimensional  analysis  based 
on  large  eddy  scale  (k,s),  it  follows  that 

[u]  =  k‘;  [£)  =  k1 '  5/e ;  [t]  =  k/e;  [£2/t]  =  k2/e 

The  (k,e)  scale  is  chosen  here  instead  of  the  small  eddy  scale  (v,e) 
based  on  the  physical  ground  that  diffusion  of  any  quantity  by  turbulent 
fluctuation  is  largely  controlled  by  large  eddy  motion.  Thus 


- P~ - - -  k2  3uiUi 

-  U.U.U.+  -(5..U.  +  6..U.)  =  C,  —  - - 

ijZ  p  jH  i  li  ke 


Here,  is  a  proportionality  coefficient.  It  should  be  remarked  here 
that  the  model  observes  the  symmetry  of  the  original  form  between  i  and 
j  as  stated  in  postulate  5.  The  dissipation  term  in  equation  (2.6)  is 
modelled  based  on  postulate  3  as 


3u .  3u .  „  9u  3u 

2  _ i  _ 1  _  2g  _ m  _ m 

dxl  *Xl  3  ijV3xfc  3x2 


e 


This  is  based  on  the  understanding  that  the  larger  the  Reynolds  number 
the  smaller  the  turbulent  eddies  are  and  that  the  smaller  these  eddies 
become  the  more  isotropic  they  will  be.  Thus,  the  dissipation  of 
turbulent  stress UTiT  by  the  small  eddies  is  mainly  in  the  isotropic 
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range.  It  should  be  noted  here  that  under  postulate  3  and  the  model 
presented,  the  dissipation  of  u^u^  can  occur  only  in  the  normal  stress 
u.u  for  i=j  and  not  the  shear  stress  when  iri  However,  when  i=i  the 

i  j  -  J 

model  term  reduces  to  the  exact  expression. 

The  pressure-strain  term  is  modelled  based  on  postulates  4  and  5  as 
[19] 


"du,  au. 


pfe  +  =  -C1  k(uiuj  -  shjV  '  C2(Pij  ' 


where  and  are  model  coefficients  determined  from  experiments  and 


_  au.  _  au.  au 

P  -  .  =  “  (u.u..  t— ^  +  u.u..  r  and  P.  =  -  u  u  t— ^ 
ij  l  1  3x.  j  1  3x/  ■  k  n  m  3x 

i  .1  m 


Further  details  of  the  modelling  of  u^uj  equation  can  be  found  in 
[18,19].  The  modelled  "uTin  equation,  thus,  has  the  form 


Du  .u . 
_ i_l 

Dt 


i_ rC  £ 

8*t‘  k  « 


2  3u  ,u 


3x 


3u  ,u . 

i— 1  +  v— 


ax. 


]  -  P. .  |5. ,e 

ij  3  ij 


C1  stVj  -  fv1  -  c2ipij  -  fvu 


(2.10) 


From  this  equation,  the  k-equation  is  obtained  by  summing  i=j  for 
i=l,2,3  and  dividing  the  result  by  two.  This  gives 


Dk 

Dt 


9  rr  k2ak 
ax£l  k  e  3x£ 


+  v 


ak 

3x£ 


] 


_  au. 

-  U  ,Urt  T - 

i  £  3x„ 
£ 


z 


(2.11) 
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It  should  be  remembered  here  that  in  equation  (2.11),  where  i=j  ,  the 


pressure-strain  term  (p/p) [ Su^/ax^  +  au^/ax^j  is  identically  equal  to 
zero  due  to  incompressibility  requirement.  Therefore,  in  equation 
(2.11),  only  the  first  terra  on  the  right  hand  side  is  modelled  and  the 
rest  of  the  equation  is  exact  as  derived  in  equation  (2.7).  It  should 
also  be  noted  that  equation  (2.11)  portrays  the  interaction  of  all 
turbulent  eddies.  The  last  term  in  equation  (2.11),  z ,  is  dominantly 
associated  with  the  small  eddies  and  is  responsible  for  dissipation  of 
turbulent  energy  that  is  produced,  first,  by  -uYu  (aU^/Sx^)  through  the 
stress  exerted  by  the  fluctuating  motion  on  the  mean  flow  motion  and 
secondly,  by  turbulent  and  viscous  diffusion  shown  in  the  first  term  on 
the  right  side.  The  diffusion  term  can  be  reasoned  to  be  more 
intimately  correlated  with  the  large  eddy  motion.  This  is  why  the  length 
and  time  scale  of  large  eddies  [1]  =  Vl1'*/z  and  [t]  =  k/e  is  adopted  in 
modelling  the  diffusion  term.  Although  two  scale  concept  is  evident  in 
the  k-equation,  there  is  no  need  to  invoke  the  second  and  small  scale 
(e,v)  in  this  equation  as  the  last  term,  e,  is  exact.  The  situation, 
however,  is  different  when  one  attempts  to  model  the  e-equation .  This  is 
considered  in  the  following  section. 


2.3.2  Modelling  of  z -equation 

The  modelling  of  e-equation  is  important  because  it  governs  the  way 
in  which  the  turbulent  kinetic  energy  is  dissipated.  As  mentioned 
earlier,  the  performance  of  the  modelled  e-equation  based  on  a  single 
turbulent  scale  of  large  eddies  is  not  as  satisfactory  as  the  other 
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modelled  equations.  First,  the  model  constant,  which  appears  in  the 


equation  for  eddy  viscosity  and  Cgl>  which  appears  in  the  e -equation 
are  found  not  to  be  constants.  Secondly,  the  prediction  of  turbulent 


flow  is  quite  sensitive  to  the  values  of  the  constants  C  and  C 

el 


In  modelling  the  e- equation,  equation  (2.8),  it  should  also  be 
remarked  that  all  eddy  motions  contribute  in  the  equation.  The 
dissipative  action  is  dominant  at  the  small  eddy  level  while  the 
convective  and  diffusive  actions  are  predominant  at  the  large  eddy 
level . 

The  scale  at  which  the  small  eddy  is  manifesting  its  dissipating 
function  in  a  given  flow,  is  intimately  related  to  the  large  scale 
structure  and  the  ratio  of  inertia  force  and  viscous  force  or  the 
Reynolds  number  as  already  discussed  in  section  2.2.  The  effect  of 
large  scale  motion  on  the  small  eddy  scale  is  transmitted  through  the 
transfer  mechanism  created  by  the  nonlinear  term  of  the  transport 
equation.  Each  term  in  the  t -equation  contributes  differently  in  a 
different  range  of  eddy  size.  Thus,  it  is  important  to  model  each  term 
in  z -equation  individually  according  to  the  eddy  size  that  characterizes 
the  physical  process  of  the  term.  Proceeding  in  this  way,  the  e- 
equation  is  modelled  below.  The  first  production  term  in  equation  (2.8) 


3U .  3u.  3u.  3u.  3u . 

l ,  Z  l  ,  i  i 


3x .  3x  .  3x .  3x.  3x. 

J  i  J  l  l 


is  modelled  as 
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This  is  because,  for  i= j ,  the  mean  strain  is  zero  from  conservation  of 
mass  for  incompressible  flow  and  for  irj,  the  quantity  in  the 
parenthesis  is  zero  from  the  isotropic  nature  of  small  eddies  at  large 
Reynolds  number  which  is  mentioned  in  postulate  3. 

The  second  production  term. 


3u.  a2u. 

1  i 


VUJt  3x.  3xrt3x. 
j  *  J 


is  also  neglected,  based  on  Lumley’s  [17]  proposal  that  the  correlation 
coefficients  between  two  quantities,  each  from  a  different  range,  are  of 


another  range  by  the  small  eddy  scale  and  consequently  the  value  of 
correlation  coefficient  is  considered  small  or  weak  compared  with  the 
other  terms  in  the  £ -equation.  For  example, 


3u . 5u7 


j  £  £ 


has  a  strong  correlation  as  the  terms  3u^/9x^,  Su^/Sx^  and  3u^/9x^  are 
in  the  same  range.  Therefore, 


So  the  second  production  term  is  dropped  from  the  £ -equation. 
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The  modelling  of  the  two  destruction  terms  is  done  based  on  postulate 

4  that  they  are  function  of  the  quantity  P^/e  and  other  transport 

variables  in  accordance  with  Lumley's  argument  [17].  Here,  P  [= 

k 

”uiuj (3Ui/3xj) ]  is  the  production  of  turbulent  kinetic  energy.  Thus 


2v 


3u .  3u  .  3u . 

i  i  l 


3x.  3x„  3x„ 
j  Z  Z 


"FuT 


2[v 


1  ]  =  fn(^’  k,  s,  v) 


3VX£ 


Lumley  assumed  that  these  two  terms  should  vanish  when  the  turbulent 
flow  approaches  equilibrium.  Thus,  for  small  deviations  from  turbulent 
equilibrium,  this  function  may  be  approximately  expanded  to  obtain 


fn(r>  =  ifni  •  h  «[£][«  -  pk] 

Here  [e/t]  is  the  dimension  needed  so  that  the  overall  dimension  of  the 
s -equation  and  that  of  the  two  destructive  terms  are  consistent,  t  is 
the  time  scale  that  characterizes  the  physical  action  for  destruction  of 
e.  Since  the  dissipation  or  destruction  of  z  physically  is  dominated  in 
the  small  eddy  range,  the  time  scale,  t=(v/s)h  based  on  Kolmogorov 
hypothesis,  is  used  in  the  present  work.  Hence,  the  modelled 
destruction  term  is 


2v 


3u .  3u.  3u, 
l  i  l 


3x.  3x.  3x. 
J  l  l 


-  2[v 


"Fu. 


]  = 


Cel(e/V)*uiuj  3^T  '  CE2(e/v)te 


3U . 

l 

j 


This  model  differs  from  the  existing  turbulence  model  [19]  in  that  the 
Kolmogorov  scale  (s,v)  is  used  for  the  scaling  of  time  [t]  instead  of 
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the  conventional  scale  based  on  (k,s)  which  leads  to  the  conventional  t- 
equation  given  in  equation  (1.10). 

The  diffusion  term  in  s-equation  is  modelled  according  to  postulate  2 
that  it  is  proportional  to  the  gradient  of  s,  or 


2v  jUi  3P_  _  _  , 1 3e  =  k2  3s 

p  3x.  3x.  st  3x ,  s  e  3x. 
J  J  J  J 


Here  the  length  and  time  scales  are  modelled  based  on  the  large  eddy  or 


l  =  k1-5/E  and  t  =  k/s 


as  shown  before.  Thus  the  modelled  s -equation  based  on  the  two-scale 
concept  and  Lumley's  suggestion  for  destruction  term  is 


Ds  _  3 _ k2 3s  3e_ 

Dt  3x^  e  e  3x^  V3x^ 


,  au.  , 

C£1(£/v)  u.u.  -  Cej(£/v)*e 


(2.12) 


Here  Cg ,  and  are  model  proportionality  coefficients.  In 

general,  they  can  be  a  function  of  fluid  or  flow  properties  such  as 
Prandtl  number  or  Reynolds  number. 


2.3.3  Modelling  of  u  0 -equation 


Finally  the  modelling  of  the  u^0- equation  [equation  (2.9)]  is  done  to 
complete  the  turbulence  closure  problem.  Using  the  same  postulates  as 
those  for  modelling  the  u_^u^  equation,  the  diffusion  term  of  this 
equation  is  modelled  as 
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u  .u.  0 

1  Z 


'  TV  ,  2  3u  .  9 

5  1I=  c 

Up  T  z  3x 


and  for  a  =  v, 


30 


3u .  3u . 0 
,  n  l  i 

au .  - —  +  v0t —  =  a; - 

1  3xi  Sx«  3x» 


The  dissipation  term  vanishes  due  to  the  assumption  of  isotropic  nature 
of  small  eddies  or  postulate  3,  i.e., 


-  (a  +  v)t —  - —  =  0 

sxt  axt 


The  pressure  strain  term  ( P — 8 )  is  modelled  according  to  Launder  [20]  as 


P  30_ 
p  3x . 


, _  au. 

-  C  -u.0  +  C  u  0-  ~ 
T1  k  i  T2  m  3x 

m 


The  frictional  term  in  "O’- equation  is  neglected  as  it  is  an  order  of 
magnitude  smaller  than  the  other  terms.  Hence  the  modelled  TO- equation 
takes  the  form 


DiO 

l 

Dt 


3u .  0 


1 _ rr  k  i  .  1  , 

3x.  ^T  e  3x  +  “a-  1 


a  3xn 

_  au. _ 

-  C  ru  .  0  +  C_,  T-^u  0 
1 1  k  i  T2  3x  m 
m 


aT  3U. 

■Vi  TT  -  Vsri 

Z  Z 


Again,  C^,,  and  C^9  are  model  coefficients. 


(2.13) 
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Summarizing  the  two-scale  second  order  turbulence  model,  we  have  the 
following  equations 


3U. 

_ 1 

3x. 


=  0 


(2.1) 


du,  aD  at.,  at..' 

p_i  =  G  .  +  — n  +  -ii 

Dt  H  i  3x.  3x.  3x. 


(2.2) 


DT 


pCDt  =  1J 


3U.  3q.  aq 

T . .  1  1  nl 


3x.  3x .  3x. 

J  i  i 


(2.3) 


Du.u. 

Dt 


.  ,2  3u,u,  3u .u . 

3  k*  i  j,  +  v_i_j 


3x.  e  3xa 
l  l 


3x  ^  Pij  35ij£ 


C.  J[u.u.  -  ffi.  .k]  -  C0[P..  -  h.  ] 
1  kl  i  j  3  ij  J  2l  ij  3  ij  kJ 


(2.10) 


Dk  _  3  k2ak  ak  _  9Ui 

Dt  -  a^fck  ?  aT-  +  '-3^1  -  “i“t  JT  '  ‘ 
III  l 


(2.11) 


§t  -  f^ICE  E  Hq  +  •  cti(E/v)iuiuj  a^  -  ce2(£/'')1' 


au. 


(2.12) 


Du.e 

1 

Dt 


a _ ,  _  k2S+  2i 

3x£  Ct  e  ax£  “ax^  ] 


aT  au- 

[uiu*  a^  +  Va^ 


au._ 

•  CT1  k“?  +  CT2 

m 


(2.13) 
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2.4  Determination  of  Turbulent  Coefficients 

The  above  11  equations  for  turbulent  quantities  have  9  dimensionless 
coefficients  or  constants,  namely,  Cfc,  C£,  Cg ,  Cgl,  Cg2>  Qj.,  CT  and 
CT2  t0  be  determined  from  experiments.  The  determination  of  the  value  of 
turbulent  coefficients  in  principle  is  similar  to  the  one  for  laminar 
flow  where  an  experiment  has  to  be  performed  to  obtain  the  values  of 
viscosity  and  thermal  diffusivity  of  the  fluid.  The  laminar  coefficients 
which  are  dimensional  such  as  kinematic  viscosity  v  and  thermal 
diffusivity  a  turned  out  to  be  dependent  on  fluid  and  thermodynamic 
variables,  temperature  T  and  pressure  P.  The  turbulent  coefficients  are 
dimensionless  and  can,  in  general,  be  functions  of  fluid  and  flow 
properties  such  as  Prandlt  number  or  Reynolds  number.  If  the  turbulent 
flow  equations  are  properly  modelled,  the  model  coefficients  should 
remain  universal  and  can  be  evaluated  once  for  all  from  the  chosen 
experiments.  Thus,  the  process  of  determining  the  constants  is  not  a 
case  of  experimental  data  fitting.  It  should  be  remarked  that  although 
these  coefficients  may  depend  on  fluid  properties  like  laminar  flow 
coefficients  v  and  a,  they  are  determined  mainly  from  experiments 
performed  in  air  and  water.  Many  investigators  consider  that  these 
coefficients  remain  the  same  for  both  fluids.  Whether  these 
coefficients  are  valid  for  turbulent  flows  in  other  fluids  such  as  oil 
or  liquid  metal  is  not  known.  The  following  subsections  highlight  the 
method  of  obtaining  these  constants. 
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2.4.1  C 


and  C 


The  coefficients  and  ^e  obtained  from  experimental  data  of 
homogeneous  shear  flow  and  turbulence  behind  a  grid  [21],  Consider  a 
uniform  flow  of  velocity  Uq  passing  a  square  grid  with  spacing  M.  Tb* 
flow  behind  the  grid  can  be  made  isotropic  by  contracting  the  area  of 
cross  section  by  a  factor  of  1.27.  The  k  and  s  equations  for  isotropic 
turbulent  flow  behind  a  grid  (figure  2.3)  are 


and 


Here  x  is  the  coordinate  along  the  flow  direction.  It  should  be 


remarked  here  that  the  diffusion  terms  of  k  and  e  equations  in  their 


exact  form  are  zero  for  isotropic  flow.  Nondimens ionalizing  these 
equations  using  the  variables 


U  M 
o 


o  .  o 


V 


the  following  equations  are  obtained,  i.e., 


and 
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as  ' 
dx' 


-C  Re  ~s 
z2 


5 


From  figure  2.4,  the  relation  between  k  and  x  is  found  to  be 
approximately  for  air  [18] 


k* 


450  ,-1 

40000  X 


between  10  <  x/M  <  200  for  the  approximately  isotropic  range  and  for  Re 

M 

3  4 

ranging  from  10  to  10  .  Substituting  this  in  the  k-equation  gives  a 
relation  between  e  and  x,  which  is 


dk'  _  450  ,-2 

dx'  40000X 


This  is  now  substituted  in  the  s -equation  to  give 


450 

40000 


2  = 


c  nr*f..450 

Cs2R  40000 


)1-5x' 


-3 


Hence,  C  is  calculated  to  be 
Z2 


r  -  gr40000  j  -j_  18.9 
s2  2(  450  5  Re  ~  7rT 


It  should  be  mentioned  that  the  flow  behind  uniform  grid  is  not  truly  an 
isotropic  flow  since  u2/v2  is  always  greater  than  one.  u2/v2  starts  with 
about  1  immediately  behind  the  screen  for  Re^  >  103  and  increases  to 
1.55  downstream  [18].  Therefore,  decay  data  for  turbulent  kinetic  energy 
k  versus  x  beyond  x/M  >  200  should  not  be  considered  as  an  isotropic 
data  and  used  in  determination  of  turbulent  coefficients.  The  value  of 
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is  next  obtained  from  homogeneous  shear  flow  as  shown  in  figure  2.5. 
For  such  a  flow,  the  k  and  z  equations  are 


n  -  3U 
0  =  -uv—  -  z 

8y 


and 

0  =  “C,!**/")*  uv|^  -Ce2(e/v)^e 


Substituting  the  k-equation  into  the  £ -equation  gives 

C£l(e/v)^E  -  C£2(s/v)^e  =  0 

Therefore, 


C 


El 


18.9  Re 


-i 


The  diffusion  terms  in  k  and  e  equations  are  assumed  to  be  approximately 

zero  here.  Strictly  speaking  they  are  nonzero.  Consequently  the 

determination  of  should  be  considered  only  an  approximate  one.  The 

coefficients  and  C£2  for  the  destruction  term  in  the  £ -equation  were 

found  to  be  function  of  Reynolds  number  based  on  a  characteristic  mean 

flow  velocity  Uq  and  a  characteristic  length  M.  The  appearance  of 

Reynolds  number  in  C  .  and  C  „  reflects  that  the  small  eddies 

zl  z2 

responsible  for  destruction  for  z  are  indeed  a  function  of  mean  Reynolds 
number.  In  other  words,  the  size  of  small  eddy  and  the  time  scale  that 
characterizes  the  destruction  of  e  changes  when  Reynolds  number  changes. 


47 


It  should  be  remarked  here  that  in  the  one-scale  turbulence  model  the 

coefficients  and  in  equation  (1.10)  are  found  to  be  independent 

of  Reynolds  number.  Their  values  are  not  universal  and  require 

modification  in  some  flow  configurations  such  as  between  plane  jet  and 

round  jet.  To  compensate  the  contribution  of  diffusion  C  in  the 

£  1 

present  study  is  taken  to  be  approximately  17.5  Re'2.  The  fact  that  C 

E  1 

and  Ce2  required  modification  weakens  the  predictability  of  the  one 
scale  turbulence  model  and  motivates  the  present  investigation  of  the 
two  scale  turbulence  model  to  improve  the  predictability  of  the  model. 


2.4.2  C  and  C 


The  constants  and  C 2  in  equation  (2.10)  are  obtained  in  a  way 
similar  to  and  C^2  [18].  Experimental  result  of  anisotropic 
turbulence  behind  a  grid  by  Uberoi  [21]  are  used.  For  such  a  flow, 
U=Uo=constant  and  V=W=0 .  By  passing  this  flow  through  a  4:1  contraction 
of  flow  cross-section  area,  the  turbulence  becomes  strongly  anisotropic. 
With  Uq=  constant,  the  exact  equation  for 'uTuT  [equation  (2.6)]  where  TT3 

""  x  - r 

>  v  =  w  becomes 


U 

o 


o  3u  3u  2p  3u 
2v~ —  - —  +  - — 

3Xfc  9X£  P  3x 


when  i— j  and  the  modelled  equation  is 


U 

o 


Is  •  ci  i(u 


t  ,-~2  2. 

-fir  -  — 


k) 


2/u2  (*104 ) 
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Figure  2.3.  Isotropic  flow  behind  a  grid 


Dividing  the  second  term  on  RHS  of  each  equation  by  the  first  term  on 
the  RHS  gives 
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2p  3u 
P  3x 


C  (— 
ClC2k 


1) 


In  order  to  accurately  determine  one  needs  to  avoid  data  of  1.5(iP/k) 
which  is  close  to  one  since  it  will  make  the  right  hand  side  of  the 
above  expression  zero.  In  other  words,  one  should  consider  the  data  in 
the  strongly  anisotropic  range,  ~ui/vr>  1,  or  between  x/M  equal  to  0  and 
40.  From  figure  2.6,  ("v2/ur)=l .  83  for  x/M=25  where ~w^=v^.  Thus 


3  _  3^*71.83 

2  k  v^ (1/1 . 83  +  2) 


0.644 


From  figure  2.7,  at  x/M  =  25, 


2P  3u 
p  3x 


2v 


3u  9u 

3x„  3x. 
I  l 


1.0 


Hence 


-  C 


3? 


l(2k 


1)  =  1.0 


and 


C1  =  2.8 

If  data  at  x/M  =  12.5  is  used,  =  2.88. 
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X/H 

Figure  2.5.  Homogeneous  shear  flow 
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Figure  2.6.  Experimental  data  of  Uberoi. 
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Next,  C?  is  obtained  from  experimental  data  of  homogeneous  shear 
flow.  Figure  2.5  shows  the  values  of  k,  u7,  v7  and  w7  for  such  a  flow. 
The  modelled  Reynolds  stress  equation  (2.10)  for  "u7  becomes 


0=0-  2uv|^ 
3y 


C1  '  fk> 


C2(2  -  fjTTvfS 


The  k-equation  (2.11)  becomes 


0=0-  uvf^  -  E 
3y 

From  figure  2.5,  at  x/H=10, 

(u1  -  |k)/k  =0.22 

Substituting  the  above  value  and' s  in  the 'll7  equation  gives 

4(1  -  C2)/3C1  =  0.22 

For  C^=2.8,  is  found  to  be  approximately  0.54.  The  commonly  used 
values  of  and  C are  2.3  and  0.4. 


2.4.3  CT1  and  CT2 

For  these  coefficients,  the  experimental  data  of  homogeneous  shear 
flow  with  a  temperature  gradient  obtained  by  Webster  [22]  is  used.  The 
modelled  equation  for  u . 0  is 
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Figure  2.7.  Experimental  data  of  Uberoi. 
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a5?, 

'  3xe  T  1  S  X 


Dt 


_  2*7'  3U  , 

•  <vt  »r  +  ".'sri 


30. _ 

-  c_,  rO  +  r  — — -u  0 

T1  k  i  T2  3x  m 
m 


(2.13) 


To  determine  ,  i  is  set  to  2.  For  the  homogeneous  shear  flow,  this 
gives 


0  =  0  -(pg  -  0)  -  cT1  fve 


and  the  k-equation  gives 


0=0-  uvjp 
9y 


Z 


Therefore, 

c  =  -  P-  (-3T/3y)  =  i  k  (3T/3y) 

T1  z  vH  uv  W  (3U/3y) 

From  experimental  data  of  Webster  as  shown  in  figure  2.8  the  magnitude 
of  uv,  v7 ,  k  and  v0  are  found  to  be  about  0.5,  1.9,  3.23.  and  0.38, 
respectively.  These  values  are  for  Richardson  number,  Ri,  of  0  as 
indicated  by  the  dashed  lines,  which  represent  the  averaged  value  of  the 
experimental  data.  Further,  the  ratio  of  the  temperature  gradient  to 
the  velocity  gradient  is  obtained  from  experimental  date  to  be  0.1. 
Substituting  these  values  in  the  above  relation,  CT  is  calculated  to  be 
3.2.  Similarly,  is  found  by  letting  i— 1  in  the  modelled  u^0 

equation.  This  gives 
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-cn  ^»*cn58|® 

Again,  from  experimental  data  of  figure  2.8,  uB  is  0.47  at  Ri  =  0. 

Hence,  is  calculated  to  be  0.5. 

2.4.4  Ck,  C£  and  CT 

The  coefficients  C^,  C£  and  C^,  are  obtained  by  computer  optimization  to 
be  0.9,  2.00  and  0.13  respectively.  Several  investigators  [18]  have 
obtained  these  constants  from  experimental  data  of  near  wall  turbulence. 
Their  values  are  not  used  in  the  present  model.  Instead,  the  modified 
values  that  give  best  results  are  used.  However,  once  the  values  are 
determined  they  are  kept  constants  for  all  calculations  in  the  present 
study. 


2.5  Concluding  remarks 

From  the  above  discussion,  the  9  turbulent  coefficients  or  constants 
are  determined  to  be  as  follows: 

Ck=0.9;  C£=2 . 00 ;  CT=0.13 

C  x=2 . 8 ;  C£l=17.5/(Re)*;  CT1=3.2 

C2=°-4;  Cg2=18.9/(Re)*;  CT2=0.5 

where  Re  is  the  Reynolds  number  based  on  the  problem  characteristic 
velocity  and  length.  These  coefficients  are  determined  from  different 
experiments.  However,  if  the  turbulence  model  is  to  have  predictability, 
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Figure  2.8.  Experimental  data  for 
homogeneous  shear  flow  (Webster) 
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the  coefficients  should  remain  the  same  in  other  turbulent  flows.  This 
will  be  examined  in  Chapter  V. 

It  should  be  remarked  that  Reynolds  number  appears  in  the  turbulent 
coefficients  for  the  two-scale  turbulence  model.  This  was  not  the  case 
for  the  one-scale  turbulence  model.  The  appearance  of  Reynolds  number  is 
expected  since,  as  discussed  in  section  2.2,  turbulent  flows  are  still 
Reynolds  number  dependent  and  also  because  small  scale  (v,e)  proposed  by 
Kolmogorov  for  dissipation  of  turbulent  kinetic  energy  contains 
kinematic  viscosity.  Physically,  this  implies  that  Reynolds  number 
changes  the  magnitude  of  the  destruction  of  z  and  hence  affect  the 
magnitude  of  z  and  range  of  eddy  size  that  is  responsible  for  the 
dissipating  the  kinetic  energy. 
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CHAPTER  III 

TWO  SCALE  K-I-CT2  TURBULENCE  MODEL  FOR 
BUOYANT  FREE  SHEAR  FLOWS 

Flow  patterns  in  nature  are  complicated  due  to  change  in  density 
caused  either  by  a  temperature  or  concentration  difference.  The  force 
produced  by  this  variation  is  called  buoyancy  force.  In  this  chapter, 
the  two-scale  turbulence  model  derived  in  the  previous  chapter  for 
incompressible  flow  is  extended  to  include  turbulent  buoyant  flows. 

3 . 1  Bous s inesq  * s  approximat ion 

and  governing  equations 

For  flows  where  the  density  gradient  is  not  large,  the  buoyant  force 
can  be  incorporated  into  the  governing  equations  by  making  the 
Boussinesq  approximation.  In  this  approximation  [18],  the  density 
variation  is  considered  significant  only  in  the  gravitational  term.  In 
other  words,  the  effect  of  density  difference  in  the  conservation  of 
mass,  in  the  time  rate  change  of  momentum  and  in  the  work  done  due  to 
density  changes  are  considered  negligible. 

The  governing  equations  for  such  flows  are 


(3.1) 
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a. 

DU" 
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3x , 
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32U . 

’i  +  y3x.3x. 


and 


(3.2) 


*  * 
DT  32T 

PsCDt  3x.3x. 

J  J 


* 

+  0 


(3.3) 


0  is  the  heat  source  due  to  dissipation  by  viscous  force.  c  or  c  and 

P 

K  are  the  specific  heat  and  thermal  conductivity.  The  superscript,  *, 

Vr  Vr 

represents  an  instantaneous  quantity.  If  P  and  p'  at  the  static  state 

jl 

(U^  =0)  are  and  p^,  then  the  momentum  equation  becomes 


0  = 


3P 

S 

3x 


+  PsSi  +  o 


(3.4) 


Subtracting  the  two  momentum  equations  (3.2  and  3.4),  we  have 


P 


s  Dt 


32U. 


•  ps)8i  +  “37-57: 


(3.5) 


With  the  Boussinesq  approximation,  the  pressure  and  density  relations 
are  given  by 

p*  *  P*  -  Ps;  (P*  -  Ps)/Ps  =  -  P(T*-  Ts)  =  -  SAT* 

* 

where  P  is  the  pressure  above  the  static  state  and  T  is  the  absolute 

s 

temperature  at  the  static  state.  3  is  the  volumetric  expansion 
coefficient  or  -l/(J(3p/3T)p  evaluated  at  T  ,  P  .  The  governing 
equations,  therefore,  become 
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and 
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DT  _  32T 

Dt  a3x.3x. 


*  * 
x  . .  3U . 

+  — ^  -i 

pc  dX  . 
J  J  J 


(3.8) 


is  the  instantaneous  viscous  stress  and  a  is  the  thermal 
diffusivity  or  K/p^c.  These  equations  are  exactly  the  same  as  the 
equations  for  non-buoyant  flows  except  that  the  momentum  equation  has  an 
additional  buoyancy  term  (JAT  g  .  Letting 

*  *  *  * 

4-  u^;  P  =  P  +  p;  x  ^  =  x..4*t..,;T  =  T  +  0 

and  taking  an  ensemble  average  of  these  equations,  the  resulting 
equations  obtained  are 
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and 
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The  term 


0) 


may  be  considered  as  frictional  heating  due  to  the  fluctuating  and  mean 

flow  motion.  This  is  normally  small  in  turbulent  buoyant  flows  and  can 

be  omitted  in  most  of  buoyant  flow  studies .  The  above  equations  have 

terms  u  u.  and  u . 0  which  need  to  be  modelled.  The  equations  for  u  u  k 
i  J  i  i  j 

e  and  u~U  are  obtained  in  the  same  way  as  those  in  chapter  II.  The  final 
form  of  these  equations  is 
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where 
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3u  .  3u  .  3U  .  3U  .  3U  ,  3u  .  3u  ,  3u  .  3u , 
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In  equations  (3.12)  and  (3.14),  each  of  the  terms  -3(g“u7&  +  *TTT)  and 

1  J  j  i 


-6g /uTU  is  called  the  buoyancy  production  and  is  a  new  source  term  in 

the  budget  of  Reynolds  stress  and  turbulent  kinetic  energy.  The 

turbulent  heat  flux  pcuTTT  now  assumes  an  additional  role,  because  it 

participates  in  the  production  terms  for  both  k  and  T2 .  In  the  uT 

i 

equation  above,  the  new  term  due  to  buoyancy  is  S$g . 0 *  which  needs  to  be 
modelled  as  it  is  an  additional  unknown.  To  model  it,  a  transport 
equation  for  0 1  is  next  derived.  This  equation  is  obtained  by 
multiplying  the  equation  for  0  by  20  and  taking  the  ensemble  average. 


This  results  in 
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The  quantity  IP  can  be  considered  as  the  intensity  of  temperature 

fluctuation.  Thus,  pc^/lP*  represents  the  fluctuating  thermal  energy.  In 

other  words,  IP  is  to  the  turbulent  heat  flux  as  k  or  u  ,u~ /2  is  to  the 

i  i 

turbulent  stress.  In  equation  (3.16),  the  rate  of  change  of  82  is 
controlled  by  turbulent  and  molecular  transport  of  IP  (  the  first  two 
terms  on  the  right  hand  side  of  the  equation  ),  the  gradient  production 
(  which  is  like  the  production  term  of  turbulent  kinetic  energy  ),  by 
molecular  dissipation  (  a  is  the  thermal  diffusivity  )  and  the 
frictional  heating  (the  last  term) . 


The  molecular  dissipation  of  temperature  fluctuation  (the  fourth  term 

on  RHS  of  equation  (3.16))  is  similar  to  the  dissipation  of  turbulent 

kinetic  energy,  e.  This  term,  a(30/3x^) (39/3x^) ,  is  denoted  by  in 

analogy  with  £,  and  represents  the  dissipation  of  the  temperature 

fluctuation  (P^  or  the  fluctuating  thermal  energy.  sQ  is  an  unknown  and, 

o 

therefore,  an  empirical  relation  or  a  transport  equation  similar  to  e- 

equation  is  needed  to  solve  it.  Here,  the  transport  equation  for  is 

o 

derived  by  differentiating  equation  (2.5)  with  respect  to  x^,  then 
multiplying  it  by  2o(30/3x^)  and  taking  the  average.  This  gives 
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It  should  be  noted  that  £g  appears  only  in  the  "Q^-equation  and  hence 
plays  a  less  important  role  than  £ -equation  in  the  determination  of  the 
mean  flow  quantities  such  as  P  and  T.  For  this  reason  some 
researchers  [23]  employed  a  simpler,  empirical  relation  to  model  the 
behaviour  of  e0  to  minimize  the  complexity  in  turbulence  modelling.  This 
will  be  discussed  in  more  detail  later.  For  buoyant  flows,  the  terms  u .  t) 
appear  in  both  u  and  k  equations  and  so  the  turbulent  momentum 
transfer  and  thermal  energy  transfer  are  now  coupled. 


All  nonbouyant  terms  are  modelled  the  same  as  before  in  chapter  II 
and  so  only  the  -modelling  of  the  buoyant  terms  is  given  in  the  next 
section. 


3.2  Turbulence  model 


In  equation  (3.12)  for  TTTuT,  the  only  term  that  requires  additional 
modelling  is  the  pressure  strain  term.  To  model  this  term,  the 
divergence  of  equation  for  fluctuating  velocity,  which  is  equation  (3.7) 
subtracted  from  equation  (3.10),  is  taken  to  give 
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Using  Green's  theorem,  pressure  is  obtained  to  be 


66 


E 

P0 


.1-  / 


4ir  vo£ 


(u.u 
£  ra 


u„u  )  3U„  3u 

- 2_SL  +  2— *  — S  +  - 

3x-3x  3x  3x.  P°£  3x„ 

1  m  m  £  £ 


3  0  , dvo 1 


if  rP  8  jl  I  LLi^ 

it  Js  p  3x  r  r  3xlp  S 
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where  the  surface  integral,  the  second  term  on  the  right  side,  is 
negligible  for  flows  far  away  from  the  solid  wall.  Thus,  when  the  volume 
of  integration  is  sufficiently  large,  the  pressure  strain  term  in 
equation  (3.12)  becomes 
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(3.20) 


The  superscript  *  denotes  the  term  involving  the  integration  .variable. 

Here  0..  0..  9  and  0..  «  correspond  respectively  to  the  first,  second 

and  third  terms  in  the  volume  integral.  The  reason  for  dividing  it  into 

3  parts  is  that  the  pressure  strain  is  caused  by  the  fluctuating  strain 

rate,  0..  the  mean  strain  rate,  0..  and  buoyancy,  0..  The 

modelling  of  terms  0..  „  and  0..  ^  is  the  same  as  before  and  is 
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where 
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The  third  term,  which  is  a  new  term  related  to  the  the  buoyancy,  is 
modelled  as  [20] 


where 


and 


The  RHS  of  this  equation  is  the  return-to-isotropy  part  due  to  buoyancy 
and  is  similar  to  0 . .  ,  which  is  the  return-to-isotropy  due  to  velocity 

fluctuation.  These  models  are  obtained  by  contracting  the  volume  to  a 
small  value  and  ensuring  that  when  the  flow  is  approximately 
incompressible  or  i=j  the  pressure  strain  is  zero.  Hence  the  modelled 


(3.21) 
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The  modelled  k-equation  can  be  directly  obtained  from  the  above  equation 
by  substituting  i=j  and  summing  to  give 


Dk  3 _  k2  3k  3k_ 

Dt  3x„  *■  k  e  3x„  v3xB 

I  l  l 


_  3U . 

uiut  3^  '  '  £ 


(3.22) 


It  should  be  remarked  here  that  except  the  first  term  on  the  right  side 
of  equation  (3.22),  which  is  modelled,  every  other  term  in  the  equation 
is  exact.  In  modelling  the  £ -equation  (3.15),  the  new  gravitational 
term  23g^vX3uV/ 3x ^ )  (38/ 3x^ )  is  set  approximately  equal  to  zero  due  to 
the  postulation  of  isotropic  dissipation  for  small  eddies.  However,  the 
destruction  terms,  the  last  two  terms  in  equation  (3.15),  are  modelled 
based  on  the  two-scale  concept  to  include  the  influence  of  buoyancy. 
This  is  based  on  Lumley's  assumption  that  the  destruction  term  is 


Tu .  §iT  Su  Pur, 
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Constant  * 
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ltJ  1  e 
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Here  the  production  of  k  comprises  of  both  production  due  to  shear 
force,  P^,  and  buoyant  force,  P^,  in  equation  (3.22)  and  t  is  the  time 
scale  of  destruction.  Thus,  based  on  the  two  scale  turbulence  model 
concept 


9u^  3u^  3u . 

2v3xT  3x7  3x„ 
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where 
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au. 

P,  =  -  u  .u  .  - — - 
k  11  3x 


and  P,  =  -  Bg/uTU 
b  11 


Here,  again  the  time  scale,  t,  is  based  on  Kolmogorov  scale  and  not  the 
convective  scale  (k,e)  for  large  eddies.  Therefore,  the  final  modelled 
t -equation  is 


Ds 

Dt 


a  rr  k2as  3e  .  _  ,  .  A -  3Ui 

_ — ](;  —  - —  4-  v; - 1  -  C  „  fs/v  u  u  - 

dxz  E  £  3xz  3xi  Cl  *  1  j  3xi 


-  C e2(e/v)^e  +  (^(e/v)^ 


The  exact  uTo- equation  is 
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(3.24) 


The  last  term  in  equation  (3.24)  is  the  frictional  contribution  to  uTS 

i 

which  is  normally  small  and  will  be  omitted  in  the  study.  The  only  term 
that  requires  additional  modelling  when  buoyancy  prevails  in  equation 
(3.24)  is  the  pressure  strain  terra  which  is  obtained  by  multiplying 
equation  (3.19)  by  30/ and  taking  an  average,  or 
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Shrinking  the  integral  to  a  small  volume,  this  term  may  be  approximately 
set  equal  to 


3zu„u 
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Hence,  the  modelled  u70 -equation  is 
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The  0 2 -equation  in  exact  form  is 
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The  diffusion  term  is  modelled  according  to  postulate  2  in  the  principle 
of  modelling  outlined  in  chapter  II  to  be 
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Letting  the  friction  term, 
9 2 -equation  is 
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It  should  be  remarked  that  the  £g  term  in  the  above  equation  is  an 

unknown  and  can  be  modelled  by  an  empirical  relation  such  as  equation 

(3 . 28a)  based  on  Launder's  suggestion  [20]  or  by  solving  a  transport 

equation.  For  most  practical  cases,  it  is  sufficient  to  use  such  a 

simple  empirical  relation  by  assuming  that  eq  is  a  function  of  e  or 

0 

V'tir  <3-28.a> 


where  C^is  a  dimensionless  constant  and  8 ^ /k  is  required  to  make  the 
dimension  consistent.  For  more  rigorous  modelling,  the  exact  equation 
for  Eg  can  be  derived  as 
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The  diffusion  term  is  modelled  as 
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The „ product ion  term  is  small  due  to  isotropic  dissipation  and 
incompressibility.  Hence, 
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are  also  small  due  to  isotropic  dissipation.  The  destruction  term  is 
modelled  based  on  Lumley's  assumption.  Thus 
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In  the  present  investigation  of  the  two-scale  turbulence  model,  the  time 
scale  in  the  above  equation  is  taken  to  be  that  of  the  destruction  term 
in  the  e-equation  namely  that  of  Kolmogorov  scale  [t]  =  (v/g)^. 
Destruction  of  gg  is  modelled  as 
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Therefore,  the  modelled  gg  equation  is 


3  V2  1  ^rr, 

D T  =  3^t(Ce  g  +  V)3^]  '  Cel(£/V)  U£93^  ‘  Ce2(e/v)2£8 


(3.30) 


Summarizing  the  complete  two-scale  turbulence  model  for  buoyant  flows, 
the  equations  for  the  turbulent  quantities  are 
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In  addition  to  incompressible  turbulent  flow  equation  for  u^u . ,  k,  e 
and  u  0 ,  two  additional  equations  for  IT2  and  z.  are  needed  in  the 

1  O 

buoyant  flows.  Also,  in  addition  to  9  turbulent  coeffecients  C,  ,  C  ,  C  , 

k’  1’  2* 

Ct 5  ^el’  *"e2’  CT’  CT1  ^T2  neec*e<*  £°r  nonbuoyant  flows,  6  more 
coefficients  are  needed  in  the  turbulent  buoyant  flow  prediction, 
namely,  Cy  C g3>  C^,  Cg,  and  C^. 
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3.3  Governing  Equations  for  turbulent 

free  shear  flows  with  buovancv 


For  two-dimensional  turbulent  free  shear  flows,  the  above  governing 
equations  can  be  simplified  considerably.  The  assumptions  made  in 
obtaining  turbulent  free  shear  flow  equations  are 

1.  Diffusion  in  the  direction  normal  (y  coordinate)  the  flow  is  much 
.  larger  than  the  diffusion  in  the  direction  parallel  (x 

coordinate)  to  the  flow. 

2.  Pressure  gradient  is  small  in  the  flow. 

3.  Laminar  shear  stress  is  much  smaller  than  the  turbulent  shear 
stress . 

4.  Boussinesq  approximation  applies. 

5.  Frictional  heating  is  negligible. 

With  these  assumptions,  the  mean  equations  (3.9),  (3.10)  and  (3.11) 
for  turbulent  Shear  flows  under  Boussinesq  approximation  are 


(3.31) 
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and 


(3.33) 
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where  T  is  the  ambient  temperature.  j=0  for  plane  flows  and  j=.l  for 
a 

axisymmetric  flows.  Equations  for  the  turbulent  quantities  reduced  from 
equations  (3.21),  (3.22),  (3.23)  and  (3.26)  with  the  x-direction  aligned 
to  the  gravitational  vector. 


U|uv  +  Viuv 
3x  3y 


1  3  r  j  kv2  3uv,  — 3U  Q  —r 

—  — [yJC.  —  - —  -  uv—  -  Bgv9 

j  3y  k  £  3y  J  3y  a 


y 

-  C 


^k  e  3y 
_ 3U 


1  kUV  +  C2UV3^  +  C3glIFB 


(3.34) 


TT3k  ,  T.ak  1  a  r  2n  kv^ak,  — au  U0 

+  3  i?ly  ck  3  -  uv3 y*  %'  c 


(3.35) 


r.3e  „3e  _  1  3  .  j  kvok.  _  .  ,  ,2,  — uu, 
o—  +  v—  -  —  —  [yJCo  —  —r]  +  CEl(e/v)  (-uv— ) 


3x  3y  3ylJ  z  z  3yJ 


— au, 

'V 


+  C£3(E/u)*gS®  .  C  (t/v)h 
a 


(3.36) 


T73uisr  ,  ,r3u0’  _  1  3  r  j„  kv*3W,  — 3T  — *3U 

UST  +  V  -  -  CT  -  — 1  •  UV3v 


j  3y  T  z  3y 


CT1  k^  +  CT2^3y  +  s  T 


8y  9y 

u  -  cT3)_ 
— 82 


(3.37) 


382  362  _  1  3  .  j_  kv^S?1.  .-rST  „  el2 

LtT  +  vs5 - 3  ^[y  ce  1  tr1  ‘  2v953  ‘  c8i  — 


(3.38) 


It  should  be  remarked  that  since  the  term  F2  appears  in  equation  (3.37), 
a  transport  equation  for  IT2  is  used.  In  equation  (3.38)  an 
approximation  for  Eg,  as  suggested  by  Launder  in  equation  (3.28a),  is 
made  so  that  an  equation  for  sg  is  not  required. 


It  is  necessary  to  solve  all  the  above  equations  to  complete  the 
prediction  of  turbulent  buoyant  shear  flows.  However,  a  considerable 
amount  of  computational  effort  would  be  required  to  solve  the  whole  set 
of  equations.  Therefore,  some  simplifying  assumptions  are  further  made 
to  reduce  the  number  of  differential  equations  to  be  solved.  One  such 
assumption  is  to  neglect  the  convective  and  diffusive  transport  terms  of 
uTu"j  and  “O  equations.  This  leads  to  the  following  approximated 
algebraic  relations  for  the  uv,  v7,  v0 ,  u0  and  02  terms 
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In  the  present  investigation,  these  algebraic  equations  are  solved  with 
two-scale  k  and  t  equations  in  differential  equation  form.  This 
simplified  turbulence  model  is  known  as  the  two-scale  k-e  model  and  is 
perhaps  most  practical  model  for  predicting  details  of  mean  motion  and 
turbulent  transport  properties. 
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3.4  Determination  of  coefficients 


Most  of  the  coefficients  in  the  above  equations  have  been  derived  in 

chapter  II  and  so  only  the  two  additional  coefficients  and  CQ 

To  8 1 

required  in  solving  buoyant  free  shear  flows  are  discussed  here. 


3.4.1  Coefficient  C 


01 


This  coefficient  is  obtained  from  experimental  data  of  temperature 
fluctuations  behind  a  grid.  Figure  3.1  shows  the  decay  of  82  measured  by 
Gibson  and  Schwarz  [24],  From  this  figure,  02  is  found  to  vary  inversely 
as  the  three-halves  power  of  distance  behind  the  grid. 

For  such  a  flow,  the  82 -equation  [equation  (3.38)]  becomes 
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Non-dimens ionalizing  this  equation  with  the  variables 
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the  following  equation  is  obtained. 
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In  section  2.4.1,  k  and  e  were  found  to  be 


k'  =  t4IStx'-1  and  e  1  - 


40000 


40000 


Further,  from  figure  3.1,  the  relation  between  02'  and  x'  is  found  to  be 
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Figure  3.1.  Decay  of  82  measured 
by  Gibson  et  al. 
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e2'  =  3. lx 


T  -1.5 


Substituting  these  in  equation  (3.45)  gives 
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This  results  in 


C 
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3 

4 


0.75 


the  commonly  used  value  of  C0  is  0.62.  This  is,  however,  obtained  from 
Launder  s  assumption  in  equation  (3.28a).  In  the  present  investigation, 
a  value  of  0.75  will  be  used. 


3.4.2  Coefficient  C, 


T3 


This  coefficient  is  generally  set  equal  to 
obtained  earlier  in  chapter  II  to  be  0.5.  Hence, 
for  CT3. 


[20]  which  has  been 
a  value  of  0.5  is  used 


3.5  Concluding  remarks 

In  section  3.3,  it  was  mentioned  that  due  to  Launder fs  argument,  the 
equation  for  Eg  was  replaced  by  an  algebraic  relation.  Hence,  the 
structure  of  the  turbulent  heat  flux  02  is  represented  by  only  the  (k,s) 
scale  only.  In  flow  situations,  where  Eg  might  be  an  important 
parameter,  the  need  for  solving  the  complete  £g-equation  [equation 
(3.30)]  would  be  important.  This  would  bring  in  the  influence  of  the 
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two-scale  concept  or  the  effect  of  viscous  dissipation  in  the  thermal 
dissipation  or  destruction  tQ  in  the  turbulence  model. 
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CHAPTER  IV 

REVIEW  OF  EXPERIMENTAL  WORK 

4,1  General  remark 

This  chapter  briefly  reviews  various  experimental  data  for  free  shear 
flows.  The  types  of  flows  considered  are  jets,  wakes,  mixing  layers, 
coaxial  jets  and  buoyant  jets.  Reliable  experimental  data  are  selected 
for  comparison  with  the  predicted  solution  obtained  from  the  proposed 
two-scale  turbulence  model.  Table  4.1  shows  the  definition  of  the 
spreading  rate,  S,  for  different  flows  which  will  be  used  for  comparison 
later.  This  rate  of  spread  is  a  gross  parameter  independent  of  the 
distance  x.  The  symbols  used  in  the  definition  of  S  are  shown  in  figure 
4.1.  In  addition  to  the  rate  of  spread,  detail  of  velocity  and  other 
profiles  are  given  in  the  following  section. 

4.2  Jets  flowing  into  stagnant  surrounding 

The  gross  parameter  of  importance  to  the  jet  flows  is  the  spreading 
rate,  S,  which  is  defined  as 


S 


=  5i 
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Table  4.1 

Definition  of  spreading  rate  S 


Flow 

s 

Jets 

dy^/dx 

Wakes 

CW  [dvdx] 

Mixing  Layers 

d(y0.5  "  y0.9)/dx 

Here,  y^  is  the  normal  distance  from  the  jet  axis  where  the  axial 
component  of  the  velocity  is  one-half  the  centerline  velocity  U  .  The 
spreading  rate  S  is  found  to  be  a  constant  when  the  flow  becomes  self¬ 
similar  in  the  far  region.  The  definition  of  self-similarity  or  self¬ 
preserving  is  that  the  profile  or  distribution  of  dependent  variables 
are  similar  from  one  station  to  another  and  become  identical  when  they 
are  made  dimensionless  by  the  local  reference  quantities.  It  should  be 
remarked  that  experimentally  it  is  found  that  although  both  profiles  of 
the  mean  flow  variables  and  turbulent  transport  quantities  become 
similar  the  former  usually  occurs  first.  It  is  also  found  that  the 
initial  condition  at  the  nozzle  exit  affects  only  the  near  jet  region. 
Therefore,  the  rate  of  spread  and  various  profiles  far  downstream  are 
the  same  regardless  of  the  inlet  conditions. 

4.2.1  Plane  jet 


Table  4.2  summarizes  the  important  flow  parameters  measured  for  plane 
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jot  several  investigators  [45-29].  ihe  range  of  experimental  data 
was  between  5  and  150  times  the  nozzle  width.  In  each  of  the  cases,  self 
preservation  started  at  a  different  location  downstream  of  the  nozzle. 
This  could  possibly  be  due  to  the  different  initial  conditions  such  as 
velocity  profiles  or  turbulent  intensity  level.  Nevertheless,  the  rate 
of  spread,  center-line  turbulent  kinetic  energy  and  the  maximum  shear 
stres?»  in  the  self  preserving  region  obtained  by  various  investigators 
are  about  the  same.  The  results  of  Bradbury  [25],  however,  showed  that 
s e  1  f - s  imi  1  ar i t\  W’as  not  reached  and  the  rate  of  spread  continued  to 


iUv-i ease  Devond  a  distance  of  >0  nozzle  widths.  One  reason  for  this  was 
that  in  Bradbury’s  experiment  the  surrounding  air  was  not  stagnant.  Rodi 
[4]  also  found  that  jets  flowing  into  moving  surroundings  are  only 
approximate 1\  self-similar .  However,  Bradbury  [25]  indicated  that  by 
reducing  the  velocity  of  the  surrounding  air  the  velocity  and  the 
turbulent  kinetic  energy  do  not  change  appreciably. 

1 igure  4.2  shows  the  velocity  profile  in  the  self-similar  region  of  a 
plane  jet  obtained  by  Bradbury  [25],  lb  skostad  [26],  Patel  [27],  Gutmark 
._S  ;  and  Kcbins  ,-9]  .  Except  for  a  small  region  near  the  edge  of  the 
jet  boundary,  there  is  close  agreement  between  all  measurements.  Hence, 

^ ^  ptov ides  a  gooc  test  for  a  turbulence  model.  Figure  4.3  gives  the 
prorile  of  turbulent  kinetic  energy  k  in  the  far  region  as  obtained  by 
various  investigators.  There  is  a  largo  amount  of  scatter  in  these 
results.  Gutmark' s  results  seem  to  be  inaccurate  since  some 
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Table  4.2 

Parameters  for  plane  jets 


Investigator 

Bradbury 

Heskestad 

Patel 

Gutmark 

Robins 

Nozzle 

size(cm) 

46  *  .95 

150  *  1.25 

80  *  .7 

50  *  .13 

-- 

Range 

x/D 

14-70 

47-155 

12-152 

10-150 

5-100 

Reynolds 

number 

30,000 

4700- 

37,000 

35,000 

30,000 

10,000- 

60,000 

Self -pres . 
x/D 

30 

65+ 

30 

120 

60 

S 

— 

0.11 

0.103 

0.102 

0.103 

Maximum 

Reynolds 

stress 

0.026 

0.021 

0.021 

0.024 

0.02 

Max  turb 
kinetic 
energy  k 

0.067 

0.07 

0.064 

0.077 

0.064 

abnormalities  were  reported  by  him  in  his  experiment  that  the  velocity 

decay  had  an  abrupt  change  at  x/D=65  and  the  dissipation  rate  was  only 

20  %  of  the  production  of  turbulent  kinetic  energy.  Heskestad’s  results 

indicated  an  increase  in  the  value  of  u  2/U  2  even  beyond  160  nozzle 

c  o 

widths.  Axso,  Bradbury  indicated  that  his  measurements  of  v2  >  u"^  seem 
physically  unlikely.  Experimental  results  of  Patel  and  Robins  are  also 
shown  for  comparison  purpose.  In  general,  their  measurements  are 
smaller  than  that  of  Gutmark,  Heskstad  and  Bradbury.  In  figure  4.4, 
Bradbury's  measurements  of  the  Reynolds  stress  are  shown  at  two 
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different  locations.  The  maximum  value  of  Reynolds  stress  is  0.026  which 
is  slightly  higher  than  that  of  other  investigators  quoted  in  table  4.2. 
A  plot  of  the  centerline  velocity  decay  along  the  jet  axis  is  given  in 
figure  4.5  as  obtained  by  Bradbury  et  al.  [30]  and  Van  der  Hegge  [31]. 
From  this  figure,  the  potential  core  length  or  zone  of  development  is 
estimated  to  be  about  6  times  the  exit  jet  width.  Bradbury’s  data  shows 
a  lower  decay  rate  than  that  of  Van  der  Hegge.  This  slight  difference 
in  the  two  results  could  be  due  to  different  inlet  and  free  stream 
conditions  as  it  is  known  that  by  changing  the  inlet  turbulent  kinetic 
energy,  the  length  of  the  core  will  vary. 

4.2.2  Round  j  et 

Table  4.3  shows  some  of  the  gross  parameters  obtained  by  Hetsroni  [32], 
Wygnanski  and  Fiedler  [33],  Rodi  [4]  and  Shearer  and  Faeth  [34]  for  a 
round  jet.  The  measurements  for  mean  quantities  by  Hetsroni  and  by 
Wygnanski  were  done  up  to  a  distance  of  x/D=35  and  40  respectively. 
Shearer  measured  the  flow  quantities  up  to  a  distance  x/D=510  and  the 
measurements  were  slightly  different  from  the  other  two  investigators. 

It  is,  therefore,  assumed  that  the  initial  condition  still  has  a 
significant  influence  on  the  measurements  at  x/D=40  and  the  flow  profile 
may  not  reach  the  self “similar  condition.  Hence,  Shearer’s  results  are 
assumed  to  be  more  reasonable.  For  comparison  purposes,  the  latter 
result  will  be  used. 

The  velocity  profile  is  shown  in  figure  4.6  where  there  is  a  small 
variation  between  the  measurements  of  Wygnanski  and  Shearer.  The 
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Figure  4.4.  Measured  Reynolds  stress  in 
a  stagnant  plane  jet  (Bradbury) 
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Figure  4.5.  Center-line  velocity 
decay  in  a  stagnant  plane  jet 
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difference  could  be  explained  by  the  lack  of  similarity  at  x/D=40.  This 
discrepency  in  the  result  is  further  amplified  in  the  measurement  of  the 
kinetic  energy  as  shown  in  figure  4.7.  Wygnanski  and  Fiedler  obtained  a 
value  of  0.1  for  the  centerline  kinetic  energy  at  a  distance  x/D=40. 
However,  Shearer's  measurements  indicate  a  value  of  0.08  at  x/D=510 
This  value  seems  to  be  more  reliable  and  realistic  because  the 
measurements  were  taken  sufficiently  far  away  where  the  jet  is  most 
likely  self  similar  and  there  is  no  influence  of  the  inlet  conditions. 
The  variation  of  Reynolds  stress  is  shown  in  figure  4.8.  Here,  too,  the 
measurements  of  Shearer  are  slightly  different  from  those  of  Wygnanski 
whose  data  is  for  x/D=60  and  70.  Some  scatter  of  data  is  obvious  from 
the  figure.  The  decay  of  centerline  velocity  for  a  round  jet  is 
presented  in  figure  4.9.  Along  with  the  data  of  Shearer  and  Wygnanski, 
the  measurements  of  Corssin  [35]  are  also  provided.  These  data  are  more 
agreeable,  though  the  measurements  of  Shearer  start  from  x/D=50. 

4.3  Plane  wake 


One  of  the  important  global  parameters  for  the  wake  flow  is  the 
spreading  rate  of  the  wake.  The  spreading  rate  for  the  plane  wake,  S,  is 
defined  as 


s.S2i 

w  dx 
o 


Here,  y^  is  the  normal  distance  from  the  symmetry  line  to  the  location 

where  the  x-component  velocity  U  is  (U  +  U_)/2.  U  and  U  are 

c  &  c  E 
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Table  4.3 


Parameters  for  round  jets 


Investigator 

Hetsroni 

Wygnanski 

Rodi 

Shearer 

Nozzle 
size (cm) 

2.5 

2.6 

1.29 

0.1194 

Range 

x/D 

15-40 

40-98 

62-75 

170-510 

Reynolds 

number 

— 

100,000 

« 

87,000 

— 

Self -pres 
x/D 

15 

70 

62 

— 

S 

0.0713 

0.086 

0.086 

— 

Maximum 

Reynolds 

stress 

— 

0.0165 

0.0186 

0.0195 

Max  turb 
kinetic 
energy  k 

0.101 

—  - 

0.078 

respectively  the  velocity  at  the  symmetry  line  and  the  free  stream  line, 
w  is  the  defect  velocity  or  (U_.  -  U  ) . 

Extensive  measurements  have  been  made,  over  several  decades,  in  the 
wakes  of  two-dimensional  bodies.  Data  is  available  both  in  the  near  wake 
wake  regions.  The  earliest  one  was  by  Chevray  and  Kovasznay 
[36],  who  measured  the  mean  velocity  and  turbulence  quantities  in  the 
near  region  of  symmetric  wake  of  a  flat  plate.  They  measured  the 

rate  and  obtained  a  value  of  0.062  but  this  value  was  still 
increasing  with  x.  Their  spreading  rate  measurement  did  not  agree  well 
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Figure  4.6,  Measured  velocity 
profiles  in  a  stagnant  round  jet 
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in  a  stagnant  round  jet 
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Figure  4.8.  Measured  Reynolds 
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with  chat  of  others  who  obtained  values  ranging  from  0.09  to  0.11.  The 
disagreement  is  due  to  the  fact  that  the  wake  was  formed  behind  a 
streamlined  body  which  may  influence  the  wake  structure.  Patel  [37] 
observed  that  it  takes  a  distance  of  3000  for  a  wake  to  become  self¬ 
similar.  Therefore,  a  change  in  the  initial  condition  at  the  trailing 
edge  before  wake  formation  could  influence  the  mean  and  turbulence 
quantities  in  the  wake  region  for  a  considerable  distance  downstream.  A 
physical  explanation  for  this  behaviour  is  that  unlike  jet  flows  the 
production  of  turbulent  kinetic  energy  is  low-  in  the  wake  flow  and 
dissipation  is  high.  Therefore,  the  initial  or  the  upstream  conditions 
for  the  wake  must  be  accurately  prescribed  if  one  hopes  for  a  meaningful 
comparison  between  the  prediction  and  measurement.  In  particular,  the 
turbulent  kinetic  energy  level  and  the  shear  stress  may  be  influenced  by 
the  initial  conditions  far  downstream. 

Comparisons  between  prediction  and  measurements  were  limited  to  the 
decay  of  center- line  velocity  deficit.  However,  recent  measurements  by 
Andreopoulos  [38],  Pot  [39]  and  Ramaprian  et  al.  [40]  provided  abundant 
experimental  data  for  comparison  with  the  prediction  based  on  turbulence 
models.  The  measurements  of  Andreopoulos  and  Ramaprian  were  done  in  the 
near  wake  region  while  those  of  Pot  were  in  the  asymptotic  region  for 
flow  past  a  flat  plate.  Hence,  their  data  provides  a  good  test  for  the 
performance  of  the  two-scale  k-e  turbulence  model  in  both  the  near  wake 
and  far  wake.  Table  4.4  summarizes  some  of  the  work  done  in  recent 
years.  Figure  4.10  shows  the  asymptotic  velocity  deficit  profile 
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obtained  from  the  asymptotic  theory  along  with  the  result  of  Pot.  There 
is  a  slight  difference  at  the  edge  of  the  wake  which  is  probablv  due  to 
the  fact  that  the  flow  is  not  fully  turbulent  in  that  region.  In  figure 
4. 11,  the  Reynolds  stress  profile  is  shown  and  compared  with  the 
asymptotic  solution.  Except  for  a  small  region  near  the  edge  of  the 
boundary  layer  and  around  y/y^  of  0.90  ,  the  solution  lies  on  the  data 
points.  The  measurement  of  the  centerline  velocity  deficit  variation 
with  x/ 9  is  shown  in  figure  4.12.  9  is  the  momentum  thickness  based  on 
the  velocity  profile  at  the  trailing  edge  of  the  flat  plate. 


Table  4.4 

Parameters  for  plane  wake 


Investigator 

Chevray 

Ramaprian 

Andreo 

Pot 

Kovaszny 

et  al. 

polous 

Body 

flat 

plate 

flat 

plate 

flat 

plate 

flat 

plate 

Range 

x/9 

0-207 

10-79 

0-43 

3-948 

Ree 

1580 

5220 

13600 

2940 

S 

0.062 

0.12 

Max  Rey 
stress 

0.05 

0.05 

Max  Turb 
kinetic 
energy 

0.07 

0.18 

1 
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4.4  Plane  mixing  layer 

Table  4.5  summarizes  the  experimental  data  measured  by  various 
investigators  [41-46],  The  spreading  rate,  S,  is  defined  as 


dx 

where  yQ1  and  yQ  g  are  respectively  the  normal  distance  from  the 
dividing  plane  to  the  location  where  the  x-coraponent  velocity  is  0.1  and 
0.9  of  (Uj.  -  Ug).  Both  Reynolds  stress  and  turbulent  kinetic  energy  are 
normalized  with  (Uj  -  U£).  It  can  be  seen  from  the  table  that  there  is 
a  large  variation  in  the  spreading  rate.  This  is  a  major  source  of 
concern  in  recent  years  [46].  However,  recent  data  by  Husain  and  Hussain 
[47]  indicates  that  an  isolated  mixing  layer  does  reach  a  unique 
asymptotic  spreading  rate. 

Nevertheless,  the  developing  region  of  a  mixing  layer  is  not  very 
well  understood.  This  is  due  to  complex  interaction  of  the  two  wall 
boundary  layers  and  the  two  shear  layers.  For  calculation  purposes,  it 
is  important  that  well  defined  initial  conditions  and  sufficient 
turbulence  measurements  be  available  to  characterize  the  main  features 
of  the  flow.  Also,  the  data  should  cover  the  complete  mixing  region.  At 
present,  no  totally  satisfactory  set  of  data  is  available.  However, 
some  of  the  measurements  of  the  velocity  are  shown  in  figure  4.13.  Most 
of  the  data  falls  on  one  curve  indicating  that  the  results  are  in  good 
agreement.  These  results  were  obtained  under  different  conditions  at 
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Figure  4.10.  Asymptotic  velocity- 
defect  profile  in  a  plane  wake 
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Figure  4.12.  Center- line  velocity 
defect  in  a  plane  wake 
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Table  4.5 


Parameters  for  plane  mixing  layers 


Investigator 

Wygnanski 

Liepmann 

Patel 

Sami 

Dimension 

(cm) 

51  *  18 

152  *  19 

76  *  43 

30  dia. 

Range 

x/D 

58 

90 

100 

450 

Max  Re 

465,000 

1,100,000 

1,800,000 

660,000 

S 

0.2 

0.16 

0.165 

0.163 

Maximum 

Reynolds 

Stress 

0.0091 

0.008 

0.01 

0.0109  . 

Max  turb 
kinetic 
energy  k 

0.035 

0.02 

0.0275 

-- 

the  start  of  the  mixing  layer.  Measurements  of  Albertson  et  al.  [43]  and 
Sunyach  et  al.  [44]  were  in  the  initial  region  of  a  plane  jet  while  Sami 
[45]  and  Bradshaw  [46]  obtained  data  in  the  initial  region  of  a  large 
round  jet  which  is  approximately  considered  to  be  two  dimensional.  In 
figure  4.14,  the  kinetic  energy  profiles  of  self  similar  mixing  layers 
are  shown.  Unlike  the  velocity  profile,  there  is  a  large  amount  of 
scatter  and,  thus,  it  is  difficult  to  say  which  data  is  more  accurate. 
Part  of  this  discrepency  is  due  to  the  variation  in  the  initial 
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4.5  Jets  flowing  into  a  parallel  moving  scream 

Unlike  the  flow  of  jets  into  stagnant  surrounding,  this  type  of  flow 
is  known  to  be  approximately  self-similar .  Due  to  the  presence  of  the 
moving  fluid  in  the  surrounding,  the  flow  has  two  distinct 
characteristic  regions.  Close  to  the  jet  or  the  near  region,  the 
centerline  velocity  is  much  larger  than  the  free  stream  velocity 
U£,i.e. 

U 

Therefore,  the  mean  strain  rate  is  high.  In  this  region,  the  flow 
properties  are  similar  to  that  of  a  stagnant  jet.  Far  downstream,  the 
jet  centerline  velocity  is  only  slightly  larger  than  the  free  stream 
velocity, i. e. 

U 

’  1  +  5 

where  6  is  small.  Hence,  the  strain  rate  is  weak  and  the  velocity 
profile  resembles  an  inverted  wake  velocity  profile.  This  region  is 
sometimes  termed  the  'wake  like  '  region. 

Due  to  the  change  in  the  flow  characteristic  from  large  strain  rate 
to  small  strain  rate,  jets  flowing  into  moving  stream  provide  a  good 
case  for  testing  a  turbulence  model. 
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4.5.1  Plane  j et 

Although  the  spread  of  turbulent  jets  issuing  into  parallel  moving 
streams  has  been  the  subject  of  a  number  of  theoretical  treatments  [48], 
reliable  experimental  data  on  these  flows  are  still  comparatively 
sparce. 

Figure  4.15  shows  the  velocity  profile  measurements  made  by  Bradbury 
et  al.  [30]  for  several  ratios  of  U_/ (U  -  U  ).  Since  the  centerline 

L.  O  L 

velocity  decreases  with  x,  these  ratios  effectively  correspond  to 

diffsrent  x  locations  in  the  flow  field.  All  the  profiles  coincide  into 

a  single  curve  indicating  that  the  flow  is  approximately  self-similar. 

In  figure  4.16,  the  centerline  velocity  decay  is  shown.  In  both  the 

figures,  the  ratio  of  the  free  stream  velocity  U^,  to  the  nozzle  velocity 

UVT  was  0.3. 

N 


4.5.2  Round  jet 


Figure  4.17  shows  the  plot  of  mean  velocity  profiles  at  three 
different  stations  as  obtained  by  Antonia  et  al.  [49].  All  velocity 
profiles  fall  into  a  single  curve,  indicating  that  the  mean  flow  is 
almost  self -similar.  The  ratio  U£/UN  for  this  data  is  0.3. 


In  figure  4.18,  the  Reynolds  shear  stress  is  shown  for  various 

locations  of  x/D  ranging  from  38  to  248.  There  is  a  considerable  amount 

of  scatter  at  y/y^  <  0.8.  However,  the  shape  of  the  data  curve  is 

similar  and  has  a  peak  at  about  y/y  =0.8. 

2 
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Figure  4.15.  Measured  velocity  profile  of 
a  plane  jet  in  moving  surrounding 


Uo/U. 


109 


H 


U.=(UM(UM-Ue))°-a 


0.6- 


0.4- 


■  BRADBURY 


0.2- 


100 


Figure  4.16.  Centerline  velocity  decay  of 
a  plane  jet  in  moving  surrounding 
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Figure  4.17.  Measured  velocity  for  a 
round  jet  in  a  moving  surrounding 
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Figure  4.18.  Measured  Reynolds 
stress  for  a  round  jet 
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4.6  Buoyant  jets 

The  turbulence  model  for  the  prediction  of  the  turbulent  buoyant  flow 
is  given  in  chapter  III.  In  order  to  verify  the  model,  some  reliable 
experimental  data  is  necessary.  The  existing  turbulence  models  do  r'-t 
predict  the  mean  and  turbulent  quantities  close  to  experimental  data 
unless  the  model  constants  are  altered.  Chen  and  Rodi  [50]  have 
collected  available  data  on  buoyant  jets  which  can  be  used  to  verify  the 
performance  of  the  two-scale  model.  Unfortunately,  experimental  data, 
especially  the  turbulence  quantities,  for  buoyant  flows  are  not 
sufficient  for  an  accurate  test  of  the  model. 

4.6.1  Plane  buoyant  jet  ^ 

Table  4.6  shows  the  plume  region  of  buoyant  plane  jets.  The  modified 
Grashoff  number,  which  is  the  product  of  the  Grashoff  number  and  the 
heat  flux,  ranges  from  3,900,000  to  about  966,000,000.  The  Grashoff 
number  is  defined  as 


Gr 


s(pa  -  P„)D’ 

- m - 


where 

Pa=ambient  fluid  density, 
PQ=fluid  density  at  nozzle  exit 
v  =kinematic  viscosity 
D  =jet  width  or  diameter 


The  rate  of  velocity  spread  which  is  defined  as 
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is  measured  by  the  various  investigators  to  vary  from  0.095  to  0.147. 
The  recommended  value  is  0.12.  The  thermal  rate  of  spread  is  defined  as 


S 


T 


dx 


where  y^,  is  the  location  where  the  temperature  is  one-half  that  of  the 
centerline  temperature.  From  the  temperature  measurements,  the  thermal 

rate  of  spread  has  been  obtained  by  most  investigators  to  be  around 
0.13. 


Table  4.6 


Gross  parameters  for  buoyant  plane  jets 


Investigator 

Rouse 

Kotsovinos 

Harris 

Anwar 

Modified 
Gr  No. 
*10s 

39 

470 

9660 

- 

Froude 

No. 

- 

1.4-5. 9 

4-193 

16-100 

(x/D) 

650 

43 

70 

50 

S 

.15/. 14 

.095 

- 

- 

Thermal 

Spread 

rate 

.13/. 14 

.12 

.  135 

.131 
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The  velocity  and  temperature  profiles  for  a  plane  jet  are  shown  in 
figures  4.19  and  4.20.  The  first  figure  is  for  a  pure  jet  measured  by 
Bradbury  and  Van  der  Hegge  while  the  second  is  for  a  pure  plume  obtained 
by  Rouse  et  al.  [51  ].  These  are  two  extreme  cases  of  a  buoyant  jet. 

The  centerline  velocity  according  to  Chen  and  Rodi  [50]  can  be  divided 
into  three  distinct  regions  in  a  buoyant  jet,  namely,  the  near  or  the 
non-buoyant  region,  the  intermediate  region  and  the  plume  region.  In 
all  the  three  regions,  the  experimental  data  lies  closely  to  the 
theoretical  lines  which  are  obtained  from  similarity  analysis.  Hence 
the  profiles  at  all  Froude  numbers  should  lie  between  those  of  pure  jet 
and  pure  plume.  Figure  4.21  shows  the  Reynolds  shear  stress  for  plane 
buoyant  jets  obtained  by  Ramaprian  et  al.  [52].  Their  measurements 
around  y/y^=l  shows  some  scatter.  In  figure  4.22,  the  turbulent  normal 
stress  distributions  measured  by  Kotsovinos  [53]  and  Bradbury  [30]  are 
shown.  Bradbury's  data  is  for  a  pure  jet  while  Kotsovinos's 
measurements  are  for  a  pure  plume.  It  can  be  seen  that  the  turbulent 
intensity  in  a  pure  plume  is  much  larger  than  that  in  a  pure  jet. 

4.6.2  Round  buoyant  j  et 

For  a  round  buoyant  jet,  Table  4.7  summarizes  the  gross  parameters 
obtained  by  different  investigators.  The  modified  Grashoff  number  varies 
from  10*  to  101#.  The  rate  of  velocity  spread  varies  from  0.084  to 
0.12.  The  recommended  value  is  0.112.  The  value  suggested  for  the 
thermal  rate  of  spread  is  0.1. 
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Figure  4.19.  Velocity  and 
Temperature  for  a  pure  plane  jet 
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Figure  4.20.  Velocity  and 
Temperature  for  a  pure  plume 
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Figure  4.21.  Measured  Reynolds 
stress  for  a  plume 
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Figure  4.22.  Measured  normal 
stress  and  k  for  a  plume 
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Table  4.7 


Gross  parameters  for  a  buoyant  round  jet 


Investigator 

Rouse 

George 

Abraham 

Modified 

Gr  No. 

9 

*10 

1.38 

13.3 

4.53 

Froude 

No . 

- 

.714 

1.82 

(x/D) 

75 

16 

26 

S 

.084 

.  112 

- 

Thermal 

Spread 

rate 

.098 

.104 

- 

In  figure  4.23  the  velocity  and  temperature  profiles  are  shown  for  a 
pure  jet  where  the  Froude  number  is  infinite.  The  data  is  obtained  by 
Rodi  for  velocity  and  Ruden  for  temperature.  On  the  other  hand,  figure 
4.24  shows  the  velocity  and  temperature  distributions  obtained  by  George 
et  al.  [54]  for  a  pure  plume  where  the  Froude  number  is  0.  The 
turbulent  kinetic  energy  in  buoyant  jets  is  shown  in  figure  4.25. 

Rodi  s  measurements  are  for  a  pure  jet  while  George's  measurements  are 
for  a  pure  plume.  Unlike  the  plane  buoyant  jet,  there  is  a  decrease  in 
the  turbulence  intensity  for  a  round  jet  due  to  the  presence  of 
buoyancy.  Hence,  more  experimental  data  is  needed  before  a  meaningful 1 
conclusion  of  the  accuracy  of  the  turbulence  model  can  be  made. 
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Figure  4.23.  Velocity  and 
Temperature  in  a  pure  round  jet 
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Figure  4.24.  Velocity  and 
Temperature  in  a  pure  round  plume 
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Figure  4.25.  Measured  k  and  normal 
stress  in  round  jets  and  plumes 
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4.7  Summary 

From  the  experimental  data  reviewed  above,  it  can  be  summarized  that 
though  the  amount  of  data  is  enough  to  adequately  test  a  model,  it  is 
still  not  complete.  In  most  cases,  there  is  no  experimental  data  of 
turbulence  quantities  at  the  inlet,  thereby  making  it  difficult  to 
compare  the  prediction  of  flow  near  the  initial  region  particularly  in 
the  near  region  of  the  wake  flow.  Also,  in  wake  and  jet  flows,  the 
turbulence  at  the  inlet  influences  the  velocity  and  kinetic  energy  in 
the  near  region.  Therefore,  small  discrepencies  between  experimental 
data  and  numerical  calculation  using  the  turbulence  model  may  not 
indicate  that  the  model  is  unsuitable.  Due  to  lack  of  initial  condition 
some  trial  and  error  or  guess  of  initial  turbulent  condition  during 
computation  is  necessary  in  order  to  examine  or  compare  the  experimental 
data  with  numerical  results  in  some  region  of  the  flow. 

In  chapter  V,  the  prediction  based  on  of  the  two-scale  turbulence 
model  for  nonbuoyant  flows  is  compared  with  the  experimental  data. 
Chapter  VI  contains  a  comparison  of  the  results  of  the  two  scale  model 
for  buoyant  flows  and  experimental  data  for  such  flows. 
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CHAPTER  V 

PREDICTION  OF  TURBULENT  NON-BUOYANT  FLOWS 

This  chapter  presents  the  results  obtained  by  the  two-scale 
turbulence  model  for  non-buoyant  flows  based  on  equations  (2.2),  (2.3), 
(2.11),  (2.12)  and  (2.13)  with  the  turbulent  constants  of  C  =0.9, 

^=2.8,  C2= 0.47,  C£=2.00,  Cg  ^17 . 5/ (Re)*,  Cg2=18 . 9/ (Re)* ,  CT=0.13, 
Cj^=3.2  and  C,p2=0 . 5.  The  predicted  results  are  compared  with  the 
experimental  data  discussed  in  chapter  IV.  Furthermore,  prediction  of 
the  one-scale  turbulence  model  is  shown  and  compared.  For  the  one-scale 
turbulence  model,  instead  of  equation  (2.12)  for  e,  equation  (1.10)  is 
used  and  instead  of  C  =0.9,  C  =2. 00,  C  =17.5/(Re)*,  C  .=18.9/ (Re)*  the 
values  used  are  0.225,  0.15,  1.435  and  1.92,  respectively. 

As  mentioned  earlier,  the  selection  of  turbulent  free  shear  flows  as 
the  first  type  of  flow  to  verify  the  predictability  of  the  turbulence 
model  is  based  on  the  following  considerations.  First,  there  is 
sufficient  data  available  for  comparison  for  both  mean  velocity  and 
turbulent  transport  properties.  Second,  the  pressure  gradient  in  free 
shear  flows  is  negligible  so  that  the  pressure  gradient  will  not  play  a 
major  role  in  determining  the  flow  field.  Therefore,  the  prediction  of 
the  free  shear  flow  field  is  most  sensitive  to  the  modelling  of  the 
Reynolds  stress,  turbulent  kinetic  energy  and  its  dissipation.  Third, 
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the  complexity  of  the  near-wall  turbulence  is  absent  in  free  shear 
flows.  Therefore,  the  error  in  the  approximate  treatment  of  near  wall 
turbulence  can  be  excluded  from  the  problem  and  the  accuracy  of  the 
turbulence  model  can  be  carefully  examined.  Fourth,  although  it  is 
secondary,  the  numerical  procedure  in  calculating  turbulent  free  shear 
flows  is  simpler  than  that  in  wall  shear  flows  or  separated  flows. 

5 . 1  Numerical  procedure 

The  equations  derived  in  chapter  II  for  free  shear  flows  are 
parabolic  in  nature  and  so  the  GENMIX  program  developed  by  Patankar  and 
Spalding  [55  ]  is  used.  The  program  has  been  modified  by  Chen  and 
Nikitopolous  [23]  and  Chen  and  Chen  [56]  to  include  the  governing 
equations  for  k,  z  and  0  . 

Briefly,  in  the  computation,  the  two  coordinates  chosen  are  the  x  and 
T  coordinates  instead  of  x  and  y  coordinates.  The  governing  equations 

are  transformed  from  the  x-y  coordinate  system  to  the  x-Y  system.  Thus 

—T 

the  governing  equations  for  U,  T,  k,  s  and  0  are  cast  in  the  same  form. 

T 

The  initial  conditions  are  specified  for  U,  T,  k,  e  and  0  at  x=0 .  These 
conditions  for  each  flow  are  given  later  in  the  individual  section 
describing  the  flow.  The  inner  boundary  conditions  are  the  symmetry 
conditions  for  jets  and  wake.  For  the  mixing  layer,  the  free  stream 
conditions  apply  at  the  inner  side.  The  outer  boundary  conditions  are 
zero  or  constant  velocity  and  no  turbulence. 
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The  solution  at  each  section  normal  to  the  mean  flow  direction  is 
obtained  by  using  an  implicit  method.  The  marching  step  Ax  at  each 
station  is  calculated  from  various  flow  parameters.  The  grid  size  Ay  is 
0.01  times  Ax.  A  total  of  40  points  are  chosen  in  the  cross -stream 
direction  which  is  verified  by  Chen  and  Nikitopolous  [23]  to  provide 
grid  independent  solution. 

The  calculations  were  performed  upto  a  maximum  distance  of  x/D=75  for 
jet  flows  and  x/0  of  600  for  wake  behind  a  flat  plate. 

5.2  Prediction  of  gross  parameters 

The  first  thing  to  be  concerned  with  is  the  prediction  of  gross 
characters  of  the  flow  field.  For  this  the  spread  parameter  is  chosen. 
When  a  model  is  not  capable  of  predicting  an  accurate  spread  rate  for 
free  shear  flows,  it  is  not  very  meaningful  to  examine  further  the 
details  of  flow  and  turbulent  structure  in  the  flow. 


Table  5.1  shows  the  spreading  rate  for  various  non-buoyant  flows. 

For  jets,  the  spread  rate  S  is  defined  as  the  slope  of  y^  in  the  flow 

direction,  where  y^  is  the  location  in  the  normal  direction  of  a  point 

where  the  U  velocity  is  one-half  the  centerline  velocity,  i.e.  U,=0.5U  . 

i  c 

For  wake,  S  is  the  spread  rate  times  the  free  stream  velocity,  U  ,  and 

Li 

divided  by  the  velocity  defect,  w  ,  or  (U-.  -  U  )  as  defined  in  figure 

o  h  c 

4.1.  In  the  case  of  the  mixing  layer,  the  spread  rate  is  obtained  in 
terms  of  the  width  of  the  mixing  layer.  The  width  is  defined  as  the 
distance  between  the  edges  of  mixing  layer  where  the  velocity  is  10%  and 
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90?i  of  the  free  stream  velocity.  From  Table  5.1,  it  is  seen  that  the 
values  of  S  predicted  by  the  one-scale  model  for  round  jet  and  plane 
wake  without  altering  the  turbulent  constants  are  signif icantlv 
different  from  the  experimental  data  while  the  two-scale  turbulence 
model  predicts  satisfactory  results  for  all  cases  calculated.  This 
demonstrates  that  the  two-scale  turbulence  model  indeed  provides  better 
prediction  than  the  one-scale  turbulence  model. 


Table  5.1 
Spreading  rate  S 


Flow 

Type 

Spread 

Parameter 

One 

scale 

Exp . 
data 

Two 

scale 

Round 

ili 

0.1189 

0.08 

0.081 

jet 

dx 

Plane 

jet 

*» 

0.1125 

0.11 

0.109 

Plane 

wake 

»i 

0.068 

0.098 

0.0975 

Mixing 

layer 

ii 

0.159 

0.16 

0.15 

5.3  Jets  flowing  into  stagnant  surrounding 
As  mentioned  earlier,  jets  flowing  into  stagnant  surrounding  become 
self -similar  far  downstream.  For  mean  quantities,  it  should  take  a 
minimum  of  40  diameters  to  establish  self-similarity  and  about  60  to  70 
diameters  for  turbulent  quantities  depending  on  the  initial  conditions 
at  the  jet  exit.  In  the  near  region,  the  flow  parameters  are  dependent 
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on  the  initial  conditions  like  the  velocity  and  turbulent  kinetic  energy 
at  the  jet  nozzle.  It  is  necessary  to  compare  experimental  data  and 
prediction  for  both  near  and  far  regions  in  order  to  verify  the  accuracy 
of  the  one-scale  and  two-scale  models.  Unfortunately,  most  of  the 
available  data  do  not  provide  complete  information  about  the  initial 
conditions  prevailing  at  the  nozzle  such  as,  velocity  profile,  turbulent 
kinetic  energy  or  dissipation.  This  could  lead  to  differences  in  the 
decay  of  centerline  velocity  and  turbulent  Reynolds  stress  or 
dissipation  function.  The  results  for  both  plane  and  round  jets  are 
discussed  below. 


5.3.1  Plane  jet 

In  the  present  calculations,  the  initial  condition  for  the  velocity 
is 


U=UNexp(-y2) 

The  k  and  e  initial  conditions  are 

k=0 . 06U^2  exp  ( -y2 ) 

E=0.09k1>5/H 

Here  UN  is  the  jet  nozzle  velocity.  The  6 %  and  9 %  levels  of  intensity 
are  taken  here  so  that  the  predicted  result  in  the  near  field  resembles 
closely  that  of  the  measured  data.  Similar  values  were  used  by  Chen  and 
Nikitopolous  [23] . 
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Figure  5.1  shows  the  comparison  of  velocity  profile  for  a  plane  jet 
in  the  region  of  self-preservation  at  x/D=75.  The  experimental  data 
shown  in  the  previous  chapter  is  represented  once  again  in  figure  5.1 
and  most  of  the  points  fall  in  one  line.  The  dashed  line  is  the 
predicted  result  of  the  one-scale  model  while  the  chain-dashed  line 
gives  the  results  of  the  two-scale  model.  The  agreement  is  excellent 
between  calculations  and  experiment,  except  near  the  outer  edge  of  the 
jet. 


In  figure  5.2,  the  turbulent  kinetic  energy  is  shown.  As  mentioned  in 
chapter  IV,  there  is  a  large  amount  of  scatter  between  various 
experimental  data,  where  the  maximum  k  may  vary  from  0.065  to  about 
0.084.  All  the  experimental  data  are  shown  once  again  along  with  the 
predictions  of  the  one-scale  and  two-scale  models.  The  two-scale  model 
predicts  turbulent  kinetic  energy  within  experimental  scatter  near  the 
centerline.  In  the  outer  edge  of  the  jet,  the  two-scale  model  predicts  a 
larger  k.  This  larger  Ttailr  is,  perhaps,  due  to  the  numerical 
diffusion  problem. 


The  Reynolds  stress  profile  of  a  self-preserving  plane  jet  at  x/D=75 
is  shown  in  figure  5.3.  There  is  a  slight  difference  in  value  between 


experimental  data,  one-scale  and  two-scale  models  away  from  y/y.=1.5. 
The  maximum  Reynolds  stress  obtained  by  various  investigators,  as  shown 


in  table  4.2,  varies  from  about  0.02  to  about  0.026.  The  prediction  of 
the  two-scale  model  shows  good  agreement  with  the  data  by  giving  a 
larger  peak  than  the  one-scale  model.  These  predictions  can  be 
considered  accurate  within  experimental  uncertainty. 
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Figure  5.1.  Velocity  profile  for  a  plane  jet 


131 


Figure  5.2.  k-profile  for  a  plane  jet 
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Figure  5.3.  Reynolds  stress  in  a  plane  jet 
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The  results  discussed  so  far  are  in  the  self-preserving  region  of  a 
plane  jet  where  x/D  >40.  Figure  5.4  shows  the  velocity  decay  in  the  near 
region  of  a  jet.  The  slight  difference  in  the  results  can  be  attributed 
to  the  difference  in  initial  conditions  at  the  nozzle.  Chen  and 
Nikitopolous  [23]  found  that  the  larger  the  turbulence  intensity  at  the 
jet  exit,  the  shorter  is  the  potential  core  length.  They  found  that  when 
the  initial  turbulent  intensity  is  taken  to  be  6%  of  the  mean  flow 
kinetic  energy,  the  predicted  potential  core  agrees  with  that  measured 
by  experiment  [57].  Therefore,  6%  value  is  used  in  the  present 
calculations.  The  experimental  data  of  Van  der  Hegge  [31]  and  Bradbury 
[30]  are  shown  along  with  the  prediction  of  the  two  models.  With  6% 
turbulence  intensity,  the  two-scale  model  predicts  the  core  length  less 
than  that  of  the  one-scale  model.  Hence,  it  is  assumed  that  a  smaller 
turbulent  intensity  might  be  necessary  for  the  two-scale  model  to  have 
better  agreement  with  experiment. 

The  spreading  rate  obtained  by  the  one-scale  and  the  two-scale  models 
are  0.112  and  0.109  respectively,  both  of  which  are  close  to  the 
recommended  value  of  0.11  obtained  experimentally. 

5.3.2  Round  jet 

One  of  the  major  improvements  of  the  two-scale  turbulence  model  is 
the  prediction  of  round  jet  flow  field.  As  shown  in  table  5.1,  the  two- 
scale  model  predicts  correctly  a  spreading  rate  of  0.081  for  the  round 
jet  while  the  one-scale  model,  without  varying  the  turbulent  constants 
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Figure  5.4.  Centerline  velocity 
decay  for  a  plane  jet 
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from  that  used  in  the  plane  jet  prediction,  gives  0.119,  a  45%  larger 

spread  rate  than  the  experimental  value.  It  should  be  remarked  that  in 

order  to  remedy  the  deficiency  of  the  one-scale  model  in  predicting  the 

round  jet  flow  field,  many  ad  hoc  proposals  [58-60]  were  put  forth.  Pope 

[58]  proposed  that  the  constant  C  0  be  modified  as 

z2 

C  =1.9  -  0.79X 
zZ 

where 


=  If-W—  1Y  2V 
^ 3r  *  3r ^  r 


Rodi  [4]  proposed  a  modification  of  the  constant  C  „  used  in  the  z- 

z2 

equation.  The  correction  is 


y  i  dU  dU  0.2 
Ce2  =  1.92  -  0.0667^(1^1  - 

c 

Several  other  modifications  have  also  been  suggested,  namely,  by  Morse 
[59]  and  McGuirk  and  Rodi  [60]. 

The  velocity  profile  for  a  round  jet  is  shown  in  figure  5.5.  The 
initial  conditions  for  U,  k  and  t  at  the  jet  exit  are  the  same  as  those 
of  a  plane  jet.  It  should  be  remarked  that  the  abscissa  is  now  y/x.  The 
prediction  of  the  one-scale  model  with  the  same  constants  as  for  plane 
jet  gives  incorrect  spread  rate  as  shown  in  figure  5.5. 
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The  experimental  data  plotted  in  figure  5.5  are  obtained  for  a  wide 
range  of  x/D.  Hetsroni's  [32]  measurements  were  upto  a  distance  of 
x/D=35  while  Wygnanski  [33]  measured  up  to  x/D=40.  These  measurements 
are  somewhat  different  from  those  of  Shearer  [34]-  et  al.  who  collected 
data  as  far  as  510  diameters  downstream.  Overall,  these  data  points  lall 
in  one  curve  within  a  certain  amount  of  scatter  indicating  that  the  mean 
flow  reaches  similarity  some  distance  downstream  of  x/D=40.  This  is 
further  confirmed  by  the  result  predicted  by  the  two-scale  model  where 
the  result  of  the  calculations  are  taken  at  x/D=75 .  The  predicted 
results  agree  very  well  with  the  experiment  throughout  the  whole  region, 
except  near  the  edge  of  the  outer  region  where  the  measurement  may  be 
affected  by  the  intermittency  between  laminar  and  turbulent  flow.  The 
fact  that  the  centerline  velocity  predicted  at  x/D=75  agrees  well  with 
the  data  shows  that  the  velocity  decay  along  the  axial  line  is 
satisfactory. 

Figure  5.6  gives  the  distribution  of  turbulent  kinetic  energy  in  a 

round  jet.  The  calculations  are  taken  at  x/D=75.  From  this  figure,  it 

is  clear  that  there  is  a  large  variation  near  the  center  of  the  jet. 

Wygnanski' s  data,  taken  at  x/D=40  only,  indicates  a  value  of  k  /U  2  of 

c  o 

0.1  which  is  higher  than  0.08  measured  by  Shearer.  The  latter  made 
measurements  up  to  x/D  of  510.  This  difference  in  the  centerline  kinetic 
energy  could  be  due  to  the  fact  that  turbulence  quantities  become  self¬ 
preserving  much  after  the  mean  quantities  become  self-similar .  The  one- 
scale  model  due  to  its  inability  to  predict  correctly  the  spreading  rate 
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Figure  5.5.  Velocity  profile  for  a  round  jet 
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predicts  turbulent  kinetic  energy  distribution  further  away  from  the 
experimental  data.  Using  the  one-scale  turbulence  model  without  the 
modified  constants,  a  maximum  value  of  0.095  is  obtained  at  the 
centerline  while  the  two-scale  model  predicts  a  value  of  0.067. 

Although  more  experimental  data  are  necessary  to  decide  the  level  of 
turbulence  intensity  at  the  jet  centerline,  it  seems  that  the  two-scale 
model  may  predict  a  smaller  turbulent  intensity  than  the  available  data. 

Figure  5.7  shows  the  Reynolds  stress  distribution  in  a  round  jet. 
There  is,  once  again,  a  difference  of  about  20%  in  the  maximum  value  of 
the  Reynolds  stress  as  obtained  by  Wygnanski  and  by  Shearer.  The  former 
measured  a  peak  value  of  0.0168  while  Shearer  got  a  value  of  0.020.  The 
prediction  of  the  one-scale  model  indicates  a  maximum  Reynolds  stress  of 
0.025  whereas  the  value  of  0.019  obtained  by  the  two-scale  model  is 
closer  to  the  experimental  data  of  Shearer.  In  general,  the  two-scale 
model  predicts  satisfactory  results. 

Figure  5.8  gives  the  centerline  velocity  decay  for  a  round  jet  using 
the  one-scale  and  two-scale  models  along  with  the  experimental  data  of 
Corssin  [35],  Wygnanski  [33]  and  Shearer  [34],  The  difference  in  the 
result  is  due  to  the  initial  condition  of  the  turbulent  kinetic  energy 
and  the  mean  velocity  profile.  Chen  and  Nikitopolous  [23]  showed  that 
the  initial  potential  core  length  is  a  strong  function  of  the  initial 
mean  velocity  profile  and  turbulent  kinetic  energy.  For  fixed  turbulent 
intensity  at  1.25%  of  the  mean  kinetic  energy,  the  core  length  is  about 
7.3  jet  diameters  with  a  flat  exit  velocity  profile  and  3.25  with  a 
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Figure  5.6.  k-profile  for  a  round  jet 
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Figure  5.7.  Reynolds  stress  for  a  round  jet 
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triangular  exit  velocity  profile.  In  the  present  calculation,  a 
turbulent  intensity  of  6,0  and  the  exit  mean  velocity  profile  given  by 

U=UNexp(-y2) 

was  used.  The  prediction  of  the  two-scale  turbulence  model  indicates  a 
core  length  of  about  5  diameters  which  agrees  with  the  measured  core 
length. 


5.4  Plane  wake 

The  initial  conditions  used  in  the  calculations  of  the  plane  wake  are 
U=UE(y/6)1/7 

k=0.008UE2Sin[1.57(l  -  y/5)’] 
s=0 . 09k1 *  5/6 

where  is  the  free  stream  velocity  and  5  is  the  boundary  layer 
thickness  at  the  beginning  of  the  wake  flow.  Most  of  the  existing 
turbulence  models  used  in  the  calculation  of  turbulent  flows  do  not 
accurately  predict  some  of  the  flow  parameters  of  turbulent  wakes,  in 
particular,  the  spread  rate  of  the  wake.  This  is,  perhaps,  due  to  the 
fact  that  the  turbulent  process  in  a  wake  involves  complex  interaction 
among  turbulent  diffusion,  production  and  dissipation  and  also  due  to 
the  fact  that  the  flow  in  the  far  wake  region  is  not  fully  turbulent.  As 
shown  in  table  5.1  ,  the  one-scale  model  underpredicts  the  spreading 

The  measured  spreading  rate,  defined  as 


rate  by  about  30%. 


142 


as 


o 


3 


Figure  5.8.  Centerline  velocity 
decay  for  a  round  jet 
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w  dx 
o 

is  0.098  while  the  one-scale  model  gives  0.068.  On  the  other  hand,  the 
two-scale  model  predicts  quite  satisfactorily  a  spread  rate  of  0.0975. 

Figure  5.9  shows  the  asymptotic  velocity  profile  in  the  far  wake 
region  where  x/ 0=600.  The  momentum  thickness,  8',  is  obtained  from  the 
velocity  profile  at  the  trailing  edge.  The  experimental  result  is  that 
[39],  taken  at  x/8=1000,  which  compares  very  well  with  the 
predicted  profile  of  the  two-scale  model  except  at  the  edge  of  the  wake. 
When  the  predicted  velocity  profile  is  plotted  against  y/y, ,  the  one- 
scale  and  two-scale  turbulence  models  differ  somewhat  near  the  edge  of 
the  wake.  However,  it  should  be  noted  that  the  y^  predicted  by  the  one- 
scale  model  is  30%  lower  than  experimental  value. 

In  figure  5.10,  the  Reynolds  stress  versus  y/y^  is  shown  for  the  far 
wake.  The  values  calculated  from  the  one-scale  model  differ  considerably 
from  that  by  the  two-scale  model  and  experiment.  Patel  et  al.  [37] 
improved  the  one-scale  turbulence  model  with  the  turbulent  constants 
altered  and  by  introducing  the  intermittency  near  the  edge  of  the  flow. 
Their  result  showed  some  improvement  'in  the  prediction.  However,  it  is 
emphasized  that  no  modification  of  turbulent  constants  were  needed  for 
the  two-scale  turbulence  model  in  predicting  the  wake  flow. 

Figure  5.11  gives  the  center-line  velocity  defect  in  the  near  and  far 
wake  region.  U£  and  Uq  are  the  free  stream  and  centerline  velocity 
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Figure  5.9.  Asymptotic  velocity 
profile  in  a  far  wake 
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Figure  5.10.  Reynolds  stress  in  a  far  wake 
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respectively.  The  results  of  both  the  one-scale  and  two-scale  models 
are  in  good  agreement  with  experimental  results. 

5.5  Plane  mixing  layer 

The  accurate  prediction  of  mixing  layer  flow  depends  heavily  on  the 
initial  conditions  since  the  mixing  flow  resembles  the  initial 
development  of  a  jet. 

Figure  5.12  shows  the  velocity  profile  in  a  mixing  layer  obtained  at 
x/D=5.  Here  D  is  the  width  of  the  jet  exit.  The  calculations  were  done 
from  x/D=0  to  x/D=5  with  the  initial  conditions  for  U  and  k  as 

U=0  for  y>0 

U=Uj  for  y<0 

k=0.01Uj2exp(-y2)  for  y<0 

where  Uj  is  the  initial  velocity  of  the  mixing  flow.  The  calculations 
are  presented  only  for  the  two-scale  model  since  one-scale  model 
predicts  satisfactorily  the  gross  properties.  The  two-scale  model  also 
predicts  velocity  distribution  which  is  in  good  agreement  with 
experimental  results.  The  kinetic  energy  profile  obtained  by  various 
experiments  [41-46]  varies  considerably  as  shown  in  figure  5.13. 
However,  the  model  generally  predicts  satisfactory  results. 
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Figure  5.11.  Centerline 
velocity  decay  in  a  wake 
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Figure  5.12.  Velocity  profile  in  a  mixing  layer 
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Figure  5.13.  Kinetic  energy  in  a  mixing  layer 
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5.6  Jets  flowing  into  moving  surrounding 

In  this  section,  jet  flows  into  moving  surroundings  are  considered 
where  the  jet  exit  velocity  is  larger  than  the  free-stream  velocity.  As 
mentioned  in  chapter  IV,  jets  flowing  into  moving  surrounding  are  c  ly 
approximately  self-similar .  The  flow  field  can  be  approximately  divided 
into  two  regions,  namely,  strong  jet  region  where  the  strain  rate  is 
large  and  a  weak  region  where  the  strain  rate  falls  from  relatively 
large  to  small  values.  This  weak  jet  region  is  an  important  test  case 
for  turbulence  models,  since  turbulence  process  in  this  region  of  weak 
strain  rate  involves  not  only  turbulent  production  and  dissipation  but 
also  significant  amount  of  turbulent  diffusion.  Therefore,  unless  the 
turbulent  transport  equations  are  properly  modelled,  the  predictability 
may  not  be  accurate.  The  initial  conditions  for  the  jet  are 

u=uNexp(-y2) 

k=0 . 06U^2  ®xp ( -y2 ) 

and  for  the  free  stream  are 

U=U^=Constant 

K=0 
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5.6.1  Plane  j et 

Experimental  data  of  Bradbury  and  Riley  [30]  shows  that  except  near 
the  nozzle  exit,  the  velocity  profiles  collapse  into  a  single  curve 
independent  of  the  ratio 


where  U  is  the  free  stream  velocity  and  U  is  the  centerline  velocitv. 

£*  O 

With  Uq  varying  in  the  axial  coordinate,  the  ratio,  X,  may  be  used  to 
denote  various  distances  downstream  and  to  indicate  when  self -similarity 
is  achieved.  Figure  5.14  shows  the  velocity  profiles  for  a  plane  jet. 
The  calculations  of  the  one-scale  and  two-scale  models  are  shown  along 
with  the  experimental  data.  Since  the  flow  far  downstream  becomes 
approximately  self -similar ,  calculations  are  shown  for  only  one  location 
x/D=75.  The  ratio,  of  jet  exit  velocity  to  free  stream  velocity 

used  in  the  calculation  was  3.3.  The  results  of  the  two  models  show 
good  agreement  with  the  experimental  data  of  Bradbury  and  Riley  [30], 

In  figure  5.15,  the  decay  of  centerline  velocity  of  the  plane  jet  is 
shown.  The  predicted  centerline  velocity  by  the  two-scale  model  gives 
slightly  slower  decay  rate.  Nevertheless,  the  calculations  tend  to 
reach  experimental  value  far  downstream. 
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Figure  5.14.  Velocity  profile  for  a 
plane  jet  in  moving  surrounding 
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Figure  5.15.  Centerline  velocity  decay  of 
a  plane  jet  in  moving  surrounding 
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5.6.2  Round  jet 


For  a  round  jet,  the  predictions  of  the  one-scale  and  two-scale 

models  are  compared  with  the  experimental  result  of  Antonia  and  Bilger 

[49]  as  shown  in  figure  5.16.  The  experimental  data  collapse  into  a 

single  curve  for  all  distances  beyond  x/D=38 .  Hence,  the  velocity 

profile  at  x/D=75  is  shown  where  the  flow  is  approximately  self -similar . 

The  ratio  of  the  nozzle  velocity  U  to  the  free  stream  velocity  U_ 

N  E 

chosen  is  3.3.  There  is  a  good  agreement  between  the  experimental 
result  and  the  two  turbulence  models. 


The  Reynolds  stress  is  shown  in  figure  5.17  for  x/D=75.  According  to 

the  experimental  data  of  Antonio  and  Bilger,  the  stress  at  x/D=248  is 

larger  than  that  at  x/D=152.  This  indicates  that  the  turbulence 

quantities  have  not  reached  self -similarity .  The  calculations  of  the 

one-scale  and  two-scale  models  are  shown  for  x/D  of  75.  Around  y/y,  of 

2 

0.8,  there  is  some  difference  between  the  two  models  and  the 
experimental  data.  The  predicted  Reynolds  stress  in  general  follows  the 
trend  of  the  experiment  but  gives  smaller  magnitude  particularly  near 
the  peak  or  y/y^=0.8.  The  cause  of  large  value  of  measured  Reynolds 
stress  probably  is  due  to  the  initial  conditions  where  in  the  experiment 
the  nozzle  of  the  round  jet  has  a  finite  thickness  while  it  is  assumed 
infinitely  thin  in  the  computation.  Furthermore,  an  increase  in  the 
initial  turbulent  kinetic  energy  for  the  calculation  can  not  only  cause 
steeper  decay  of  centerline  velocity  than  that  shown  in  figure  5.15  but 
also  increase  the  predicted  Reynolds  stress  in  figure  5.17. 


(U-Ue)/(U„-U£) 
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Figure  5.16.  Velocity  profile  for  a 
round  jet  in  a  moving  surrounding 
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Figure  5.17.  Reynolds  stress  for  a  round  jet 
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5.7  Sensitivity  of  the  Coefficients 

This  section  gives  a  brief  discussion  of  the  sensitivity  of  the 

coefficients  and  in  the  £ -equation.  From  homogeneous  shear  flow 

and  flow  behind  a  grid  turbulence,  these  coefficients  were  found  to  be 

18.90/(Re)^  for  the  two-scale  second  order  closure  model  of  u. u . *  k  and 

i  J 

e.  However,  the  coefficient  C£l  was  modified  to  be  17.50/(Re)^  for  the 
two-scale  k-e  model.  This  modification  was  made  because,  in  the 
calibration  of  and  C^,  the  diffusion  terms  in  k  and  £  equations  for 
both  homogeneous  and  grid  turbulence  flows  were  neglected  which  is  not 
the  case.  The  values  C  =17.5/(Re)^  and  C  ^=18 . 9/ (Re)^  were  obtained  by 
solving  the  plane  jet  flow  where  the  turbulent  diffusion  term  is 
increased  in  the  calculation. 

It  is  known  that  flow  prediction  based  on  the  one-scale  turbulence 

model  is  very  sensitive  to  the  coefficient  which  has  a  value  between 

1.90  and  1.92.  Any  value  outside  this  range  may  cause  the  prediction  to 

change  significantly.  On  the  other  hand,  for  computation  based  on  the 

two-scale  turbulence  model,  C£2  may  be  changed  from  11.90/(Re)^  to 

18 .90/ (Re)^  and  from  10.50/(Re)^  to  17.50/(Re)^,  the  prediction  is 

quite  stable  and  satisfactory  as  long  as  the  same  difference  of 

1.4/(Re)^  between  C  and  C  „  is  kept. 

si  e2  * 

Table  5.2  shows  the  spreading  rate  for  a  plane  jet  for  various  values 
of  Ctl  and  C  These  calculations  were  done  for  Reynolds  number 
ranging  from  12,000  to  120,000. 
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Table  5.2 


Sensitivity  of  spreading  rate  on 
the  coefficients 


Re 

Ctl/(Re)4 

C£2/(Re)4 

S 

12,000 

17.50 

18.90 

0.1017 

M 

15.50 

16.90 

0.1146 

tt 

10.50 

11.90 

0.1200 

24,000 

17.50 

18.90 

0.1000 

120,000 

17.50 

18.90 

0.1009 

The  spreading  rate  varies  from  0.1  to  0.12  for  a  change  in  Cg2  from 
4  4 

10. 50/ (Re)  to  18. 90/ (Re) 2  at  a  Reynolds  number  of  12,000.  Hence  it  can 
be  concluded  that  there  is  a  slight  change  in  S  even  with  an  appreciable 
change  in  the  coefficients.  This  difference  can  be  attributed  more  to 
the  GENMIX  program  than  to  the  physical  phenomenon. 

The  effect  of  the  Reynolds  number  is  also  shown  in  table  5.2  for  a 

plane  jet  flow.  The  values  of  C  ,  and  C  „  used  are  17.50  and  18.90 

si  1 2 

respectively.  The  Reynolds  number  is  changed  from  10^  to  10^.  The 
change  in  the  spreading  rate  is  again  very  small.  Thus,  it  can  be  said 
that  a  change  in  Reynolds  number  does  not  affect  the  overall  structure 
of  the  jet.  It  should  be  remarked,  however,  that  by  changing  the 
coefficients  and  keeping  Reynolds  number  fixed  is  effectively  the  same 
as  keeping  the  coefficients  fixed  and  changing  Reynolds  number.  Since 
there  is  no  set  pattern  in  the  value  of  S,  the  difference  is  due  to  the 
numerical  problem. 
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To  further  study  the  effect  of  Reynolds  number,  the  kinetic  energy  is 
calculated  for  various  Reynolds  numbers.  Figure  5.18  showsthat  there  is 
very  little  difference  in  the  kinetic  energy  profile  at  different 
Reynolds  numbers.  The  difference  shown  in  this  figure  may  be  due  to  the 
numerical  diffusion  in  the  program  which  calculates  the  flow  using 
dimensional  quantities. 
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Figure  5.18.  Reynolds  number 
effect  on  plane  jet 
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CHAPTER  VI 

PREDICTION  OF  TURBULENT  BUOYANT  JETS 

This  chapter  shows  the  predictions  of  the  two-scale  turbulence  model 
for  buoyant  jets  based  on  equations  (3.31),  (3.32),  (3.33),  (3.35)  and 
(3.36)  and  that  of  the  one-scale  model  based  on  equations  (3.31), 

(3.32),  (3.33),  (3.35)  and  (1.10).  As  mentioned  in  chapter  IV,  the 
amount  of  experimental  data  available  for  turbulent  transport  quantities 
m  buoyant  jets  is  scarce  and  insufficient  to  test  the  accuracy  of  the 
turbulence  model.  Hence,  most  of  the  comparison  between  one-scale 
model,  two-scale  model  and  experimental  data  is  confined  to  mean  flow 
quantities.  The  velocity  and  temperature  decay  along  the  jet  axis  are 
shown  for  various  Froude  numbers.  In  the  present  study,  Froude  number 
is  defined  as 


g(pQ  -  Pa)D 

where 

PQ=fluid  density  at  jet  exit 
Pa=ambient  fluid  density 
UQ=jet  exit  velocity 
D  =diameter  or  width  of  jet 
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Some  of  the  measurements  of  normal  stress  u2  and  the  kinetic  energy  are 
also  presented,  whenever  available. 

6.1  Bouvant  plane  jet 

The  exit  and  initial  conditions  for  plane  buoyant  jets  are  the  same 
as  that  of  nonbuoyant  jets,  namely 

U=UN.exp(-y2) 

k=0.06UN2exp(-y2) 

e=0 . 09k1 ‘ 5/H 

For  the  temperature,  T  and  fluctuating  temperature,  02,  the  jet  exit 
conditions  were  set  as 

T  -  V(To  - 

0*=O.O6(T  -  T  )Iexp(-y1) 

o  a 

The  calculation  procedure  is  carried  out  similar  to  that  for  the 
nonbuoyant  jet  except  that  additional  equations  for  T  and  IT2  are 
inc luded . 

The  most  significant  characteristic  to  be  predicted  by  the  turbulence 
model  is  the  temperature  or  velocity  rate  of  spread  for  buoyant  jets. 
These  parameters  are  defined  as 


S 


U 
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and 


S 


T 


dyiT 


The  recommended  experimental  value  of  the  spread  rate  of  velocity  is 
0.12  while  that  for  temperature  is  0.13.  In  the  present  study,  the 
spread  rates  predicted  by  the  two-scale  turbulence  model  are  0.11  for 
velocity  and  0.135  for  temperature.  The  one-scale  turbulence  model 
predicted  values  of  0.11  and  0.116  for  velocity  and  temperature 
spreading  rates.  In  comparison,  the  two-scale  turbulence  model  seems  to 
give  better  prediction. 

The  centerline  velocity  has  been  found  to  be  a  function  of  the 
distance  x.  The  dimensionless  grouping  for  velocity 


U 


c 


pl/3 


-1/3 


versus  distance 


(|)f-2/3(!o  '1/3 
D  pa 

was  first  derived  by  Chen  and  Rodi  [50].  In  these  dimensionless  plots, 
Chen  and  Rodi  showed  that  all  centerline  velocity  decay  of  turbulent 
buoyant  jets  can  be  collapsed  into  a  single  curve.  The  buoyant  jet  flow 
can  be  divided  into  three  different  regions,  namely,  nonbuoyant , 


intermediate  and  plume  region.  From  experimental  data  the  following 
correlations  were  obtained.  For  the  nonbuoyant  region, 
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uc  .  ,  'o''2,*'1'1 

<r  -  (d> 

N  ra 


for 


0  <  ^F'2/3(— )  7  <  0.5 

D  pa 


In  the  intermediate  region,  the  relation  is 


jr  =  2.8sr-1/3(^):>/12(|: 

N  .a 


for 


0.5  £  ^F”2/3(— )  7  S  5 

D  pa 


and  in  the  plume  region, 


UN  pa 


for 


3  <  ^2/3£>"1/3 

D  pa 


-1/4 
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where  F  is  rhe  densimetric  Froude  number.  In  this  investigation,  a 
couple  of  cases  were  selected  for  calculation.  The  Froude  numbers  used 
were  6  and  24  which  lie  between  the  extreme  cases  of  F=0  for  a  pure 
plume  and  F=«  for  a  pure  jet. 

Figure  6.1  presents  the  velocity  profiles  of  the  one-scale  and  the 
two-scale  turbulence  models  in  a  plane  buoyant  jet  for  Froude  numbers  of 
6.0  and  24.  The  experimental  data  of  Bradbury  [30]  for  a  plane  jet  and 
of  Rouse  [52]  for  a  pure  plume  are  included  for  comparison.  The 
calculated  results  were  obtained  for  x/D=40.  There  is  a  good  agreement 
between  experiment  and  calculations. 

Figure  6.2  gives  the  turbulent  quantities,  namely  k  and  u2  as 
measured  by  Kotsovinos  [53]  and  Bradbury  [25].  The  kinetic  energy  at 
the  centerline  for  a  pure  jet  is  0.064  whereas  the  normal  stress  for  a 
pure  plume  is  0.14.  Increase  in  the  kinetic  energy  or  normal  stress  for 

the  plume  can  be  explained  by  the  existence  of  positive  buoyant  force 

«  - 

which  promotes  the  generation  of  turbulence.  In  other  words,  g(u9/T  )  in 

a 

equation  (3.35)  is  positive.  The  calculations  of  the  two-scale  model 
and  one-scale  model  are  plotted  for  Froude  numbers  of  6.0  and  24.  These 
profiles  fall  within  the  extreme  cases  of  a  pure  jet  and  a  pure  plume. 
Thus  the  prediction  can  be  considered  to  be  satisfactory. 

In  figure  6.3,  the  temperature  profile  is  shown  for  Froude  numbers  of 
6.0  and  24  as  obtained  by  the  one-scale  and  two-scale  turbulence  models. 
These  results  are  compared  with  the  experimental  data  of  Van  der  Hegge 
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[31]  for  a  pure  jet  and  that  of  Rouse  et  al.  [51]  for  a  pure  plume. 

Again  the  predicted  results  for  buoyant  jets  fall  within  the  two 
envelops  of  pure  jet  and  pure  plume  as  one  would  expect. 

The  Reynolds  stress  for  a  plane  buoyant  jet  is  shown  in  figure  6./  , 
The  one-scale  and  the  two-scale  model  results  are  shown  for  comparison. 
The  experimental  result  of  Ramaprian  and  Chandrasekhara  [52]  is  shown 
for  Froude  number  of  2.4. 

6.2  Buoyant  round  iet 

The  predicted  spread  rates  for  velocity  and  temperature  in  buoyant 
round  jets  using  the  two-scale  turbulence  model  are  0.1  and  0.115 
respectively.  The  one-scale  model  gives  these  parameters  as  0.12  and 
0.11.  The  experimental  values  of  the  velocity  and  thermal  spreading 
rates  are  0.112  and  0.1  indicating  that  the  two  models  predict 
satisfactorily  the  gross  parameters. 

Figure  6.5  shows  the  velocity  profile  of  a  buoyant  round  jet.  The 
two-scale  model  calculations  are  shown  for  the  two  cases  of  F=6  and  F=24 
while  the  one-scale  model  prediction  is  shown  for  F=6.  The  experimental 
data  of  Rodi  [4]  and  George  et  al.  [54]  are  also  shown  for  pure  jet  and 
pure  plume  respectively.  There  is  a  good  agreement  between  the 
predicted  result  and  experimental  data. 

Figure  6.6  gives  the  temperature  profile  of  a  buoyant  round  jet.  The 
pure  jet  profile  is  that  of  Rodi  while  the  pure  plume  temperature 
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Figure  6.1.  Velocity  profiles  for 
a  plane  buoyant  jet 
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Figure  6.2.  Turbulent  kinetic 
energy  and  normal  stress 
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Figure  6.3.  Temperature  profile 
for  a  buoyant  jet 
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Figure  6.4.  Reynolds  stress  for  a  plane  jet 
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profile  is  that  of  George  [54] .  The  calculations  of  the  two-scale  model 
and  the  one-scale  model  show  temperature  profiles  closer  to  the  pure 
plume  than  pure  jet. 

In  figure  6.7,  the  kinetic  energy  and  normal  stress  of  the  buoyant 
jet  is  shown.  The  pure  plume  data  is  that  of  George  and  the  pure  jet 
data  is  that  of  Rodi.  The  one-scale  and  two-scale  results  are  shown  for 
comparison.  As  in  the  case  of  the  nonbuoyant  jets  the  prediction  of 
kinetic  energy  by  the  one-scale  model  is  higher  than  that  by  the  two- 
scale  model. 


6.3  Concluding  remarks 

Though  the  comparison  between  prediction  and  experiment  was  not 
extensive,  it  can  be  said  that  the  results  of  the  two-scale  turbulence 
model  are  satisfactory.  More  rigorous  comparisons  are  necessary  which 
can  be  done  only  with  better  and  complete  set  of  experimental  data. 
Further,  the  effect  of  reducing  the  partial  differential  equations  for 
8 2  and  into  algebraic  equations  need  to  be  studied.  Since  the 

equation  for  was  replaced  by  a  relation  suggested  by  Launder  in 
equation  (3.28a),  the  influence  of  the  small  time  scales  on  the  buoyant 
jets  could  not  be  studied. 
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Figure  6.5.  Velocity  profiles  in 
a  buoyant  round  jet 
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Figure  6.6.  Temperature  profile 
in  a  buoyant  jet 
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Figure  6.7.  Turbulent  kinetic 
energy  and  normal  stress 
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CHAPTER  VII 

PREDICTION  OF  RECIRCULATING  FLOWS 

In  the  present  investigation  the  emphasis  is  placed  on  prediction  of 
turbulent  free  shear  flows.  This  is  because  the  mean  velocity  field  of 
turbulent  free  shear  flows  is  determined  by  the  Reynolds  stress  and  not 
by  the  pressure  force  which  is  likely  to  be  a  dominant  force  in  more 
complex  separate  flows.  Therefore,  any  turbulence  model  for  Reynolds 
stress,  turbulent  kinetic  energy  or  dissipation  of  turbulent  kinetic 
energy  can  be  better  tested  and  verified  by  its  capability  of  prediction 
without  the  action  of  the  pressure  force.  In  addition,  prediction  of 
turbulent  free  shear  flow  does  not  have  the  complexity  of  the  near  wall 
turbulence. 

In  the  previous  chapters,  the  two-scale  turbulence  model  has  been 

tested  and  verified  for  its  prediction  capability  in  turbulent  free 

shear  flows.  In  this  chapter,  the  two-scale  turbulence  model  is  used  to 

predict  turbulent  separated  flows.  Results  are  shown  and  compared  with 

those  obtained  by  the  one-scale  k-t  model  for  two  different  flows, 

* 

namely,  flow  past  a  backward  facing  step  and  flow  through  an  obstacle. 

7 . 1  Flow  past  a  backward  facing  step 

The  first  case  of  recirculating  flow  chosen  is  the  flow  past  a 
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backward  facing  step.  This  is  chosen  because  sufficient  experimental 
data  is  available  including  the  mean  velocity  profiles,  separation 
length  and  some  turbulent  quantities.  Stevenson  et  al.  [60]  used  laser 
velocimetery  technique  to  measure  mean  and  turbulent  quantities.  Eaton 
and  Johnston  [61]  reviewed  several  other  measurements  of  backward  facing 
step  for  the  1980  Stanford  Conference  meeting. 

The  flow  domain  which  was  measured  by  Stevenson  et  al.  [60]  and  a 
portion  of  the  grid  distribution  near  the  step  are  shown  in  figure  7.1. 
The  same  domain  is  used  in  the  calculation  using  the  one-scale  and  two- 
scale  turbulence  models.  The  height  of  the  step  is  H  and  the  distance 
from  the  step  to  the  upper  wall  is  also  H.  In  the  calculation,  the 
boundary  or  entry  conditions  were  set  at  a  distance  of  0.02H  upstream  of 
the  step  where  experimental  data  is  available.  The  exit  conditions  were 
set  at  12  step  heights  downstream  with  the  boundary  conditions  also  from 
experiment.  Near  the  wall  and  the  step,  the  grid  distribution  is 
nonuniform  with  more  nodes  in  the  vicinity  of  the  walls  where  the 
gradient  of  dependent  quantities  is  steep.  The  smallest  grid  size  is 
about  0.02H  while  the  largest  spacing  is  2. OH.  The  grid  system  for  the 
computational  domain  has  30  by  17  nodes. 

The  governing  equations  for  turbulent  flow  were  derived  in  chapter 
II.  These  equations  are  written  here  once  again.  The  mean  flow  equations 


are 
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and 
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(7.1) 


(7.2) 


To  complete  the  closure  problem,  the  Reynolds  stress  are  written  using 
the  eddy  viscosity  hypothesis  to  be 


-u  .u . 
i  J 


\i  t 


§k6.  . 
3  ij 


(7.3) 


The  turbulent  quantities  k  and  e  are  additional  unknowns.  The  k  equation 
is  modelled  as 


Dk  _  3  k2  3k  3k  .  -  3Ui 

Dt  -  3^[Ck  e  3xJ  V3x^]  -  uiu*  JTt  -  £ 


(7.4) 


The  s -equation  based  on  one-scale  model  is 


De  _  1 _ rr  k2lL_i  r  £ -  r  *2 

Dt  3x.[  t  e  3x.J  El  k  UiUj  3x .  '  Ce2  k 


(7.5) 


with  and  as  0*15,  1.435  and  1.92.  Based  on  the  two-scale 

concept  the  £ -equation  is 
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where  Cg  ,  Cgl  and  Cg2  are  2.00,  17. 50/ (Re)  ‘  and  18.90/(Re)^. 

The  numerical  method  used  for  solving  the  differential  equations  is 
the  Finite-analytic  method  developed  by  Chen  et  al.  [62]  [63].  The 
calculations  were  carried  by  the  computer  program  (FANS-1  Finite- 
Analytic  solution  of  Navier-Stokes  equations)  developed  by 
Sheikholeslami  [64].  This  program  incorporates  the  SIMPLER  algorithm 
suggested  by  Patankar  [55].  The  wall  function  [64]  is  used  for  the  near 
wall  velocity  and  turbulent  conditions.  The  velocity  profile  at  the 
inlet  and  outlet  were  specified  from  the  data  of  Stevenson  et  al.  [60]. 
The  kinetic  energy  profile  at  the  inlet  is  about  3%  of  the  mean  velocity 
squared.  The  dissipation  function  t  at  the  inlet  is  calculated  from  the 
turbulent  kinetic  energy  using  the  relation 


l  •  s 

*  =  0.09§ 

At  the  outlet,  both  the  turbulent  kinetic  energy  and  its  dissipation 
function  are  assumed  to  be  fully  developed,  i.e. 


3k 

3x 


=  0 


and 


3e 

3x 


0 


The  Reynolds  number  based  on  the  mean  inlet  velocity  and  step  height  H 
is  50,000. 

In  figure  7.2,  the  contour  plot  of  the  streamlines  predicted  by  the 
one-scale  and  the  two-scale  turbulence  model  are  shown.  The  time  step 


179 


BE-  . 5000E+5 

IMAX-  30 

JMAX-  17 

Y  1 

X 


Figure  7.1.  Geometry  of  flow 
and  grid  distribution 
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used  in  the  calculation  was  0.2  and  the  convergence  criterion  between 
the  solutions  of  two  consecutive  time  steps  was  set  at  0.0001  for 
velocities  and  0.001  for  pressure.  The  computational  time  on  Prime  850 
computer  for  both  one-scale  and  two-scale  model  prediction  was  60 
minutes  of  CPU  time.  From  figure  7.2,  it  can  be  seen  that  the  length  of 
separation  using  the  two-scale  model  is  about  7  times  the  step  height 
whereas  the  one-scale  model  predicts  a  separation  length  of  about  5  step 
heights.  Since  experimental  data  shows  a  separation  length  of  about  7 
step  heights,  the  one-scale  seems  to  underpredict  the  reattachment 
length.  The  reason  for  this  is  that  the  one-scale  turbulence  model  in 
general  predicts  a  larger  turbulent  kinetic  energy  and  hence  larger 
turbulent  eddy  viscosity  for  flow  after  a  step.  This  causes  greater 
mixing  or  momentum  transfer  resulting  in  a  smaller  separation  zone. 

Figure  7.3  gives  horizontal  velocity  profiles  predicted  by  the  one- 
scale  and  the  two-scale  turbulence  models  at  x/H=4.1  where  the  flow  is 
separated  and  has  a  region  of  reverse  flow.  It  can  be  seen  that  the 
two-scale  model  predicts  a  fuller  velocity  profile  in  the  mid-channel 
and  the  flow  is  separated  near  the  lower  wall  which  is  closer  to 
measured  values.  The  one-scale  model  predicts  much  smaller  reverse  flow 
than  the  experimental  results  of  Stevenson  et  al.  [60].  In  figure  7.4, 
the  velocity  profiles  at  x/H=7.1  are  given.  The  comparison  again  shows 
that  a  fuller  velocity  profile  is  predicted  by  the  two-scale  modei  in 
the  mid-channel  seems  to  have  better  agreement  with  experiment. 
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Figure  7.2.  Streamlines  using  one-scale 
and  two-scale  turbulence  models 
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Figure  7.3.  Horizontal  velocity 
profile  at  x/H=4.1 
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Figure  7.4-.  Horizontal  velocity 
profile  at  x/H=7.1 


184 


The  turbulent  kinetic  energy  profiles  at  x/H=7 . 1  and  9.1  are  shown  in 
figures  7.5  and  7.6.  The  results  of  the  one  and  two-scale  models  are 
presented  along  with  the  experimental  results .  Comparison  with  measured 
data  reveals  that  both  models  predict  correctly  the  general  trend  of 
distribution  of  turbulent  kinetic  energy  and  is  only  in  fair  agreement 
with  experimental  data.  It  should  be  remarked  that  the  peak  turbulent 
kinetic  energy  predicted  by  the  two-scale  model  is  slightly  smaller  in 
magnitude  than  that  predicted  by  the  one-scale  model. 


Figure  7.7  presents  the  shear  stress  at  x/H=9 . 1  for  the  two  models. 
Both  models  predict  correctly  the  location  of  the  maximum  shear  stress 
which  is  found  from  the  experimental  data  to  be  at  y/H^O^S  from  the 
bottom  wall.  Both  models  predict  larger  negative  stress  -*uTuT.  with  the 
two-scale  model  predicting  a  better  overall  trend. 

7.2  Flow  past  an  obstacle 


The  second  case  considered  in  this  chapter  is  the  flow  past  an 
extended  rectangular  plate  in  a  two-dimensional  channel.  Figure  7.8 
shows  the  stretched  geometry  of  the  flow  domain  along  with  the  grid 
distribution  for  flow  past  a  rectangular  plate  with  a  height  of  H  and  a 
thickness  also  of  H.  The  flow  domain  has  41  by  14  nodes.  The  same 
numerical  procedure  and  the  method  of  solution  is  used  as  that  for  the 
backward  facing  step.  The  calculation  is  performed  from  x/H=-ll  to 
x/H=45 .  The  velocity  profiles  at  the  inlet  and  outlet  are  the  one- 
seventh  power  law.  The  kinetic  energy  at  the  inlet  is  specified  using 
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Figure  7.5.  Turbulent  kinetic 
energy  at  x/H=7 . 1 
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Figure  7.6.  Turbulent  kinetic 
energy  at  x/H=9 . 1 
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Figure  7.7.  Turbulent  shear 
stress  profile  at  x/H=9 . 1 


the  experimental  data  of  Durst  et  al.  [65].  The  dissipation  at  che 
inlet  is  specified  using  the  relation 
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At  the  outlet,  the  fully  developed  profiles  for  k  and  z  are  assumed.  A 
portion  of  the  streamline  contour  near  the  obstacle  predicted  by  the 
one-scale  and  the  two-scale  model  for  Reynolds  number  of  17,000  based  on 
H  and  mean  velocity  U  is  shown  in  figure  7.9.  The  two-scale  model 
predicts  a  reattachment  zone  of  10H.  The  one-scale  model,  however, 
predicts  a  length  of  4H.  The  experimental  results  of  Durst  et  al.  [65] 
for  Re=17,000  show  a  separation  length  of  about  7H.  This  is  in  closer 
agreement  with  the  10H  predicted  by  the  two-scale  model. 

Figures  7.10  and  7.11  show  the  velocity  profiles  predicted  by  the  two 
models  at  x/H=4.1  and  x/H=7.1  respectively.  The  two-scale  turbulence 
model  has  a  general  tendency  to  predict  a  flatter  velocity  profile  in 
the  separation  zone  as  shown  in  figure  7.10  which  is  in  better  agreement 
with  the  experimental  data.  Figure  7.11  shows  that  the  two-scale  model 
predicts  separation  as  indicated  by  the  experiment  where  the  one-scale 
predicts  no  separation.  In  figures  7.12  and  7.13,  the  kinetic  energy 
profiles  are  shown  at  the  locations  x/H=4.1  and  x/H=7.1.  Again,  the  two- 
scale  turbulence  model  predicts  a  smaller  kinetic  energy  profile. 

7.3  Concluding  remarks 


The  results  presented  above  are  two  cases  of  recirculating  flows 
predicted  by  the  turbulence  models.  Comparison  with  experiment  shows 
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Figure  7.8.  Flow  geometry  and  grid  distribution 
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Figure  7.9.  Streamlines  using  both  one- 
scale  and  two-scale  turbulence  model 


191 


i  . 

*  / 

RE-  . 1700E+5 

.  73 

.  s 

f-  \ 

IMAX-  41 

X 

\ 

JMAX—  14 

Y 

X. 

LQC: 

.  23 

% 

I-  28 

* 

X-  7 . B7 

)\ 

* 

0 .  1 

ISTEP-  lOO 

-1 

L  * 

% 

9 

O 

<  2 

* 

TWO-SCALE 

H OR . VEL 


-1  . 


1  . 


3  . 


RE- 

1700E+5 

IMAX  — 

A  1 

JMAX- 

1 A 

LOC: 

I- 

ZB 

X- 

7 . 07 

ISTEP— 

lOO 

TWO— SCALE 


HOR . VEL 


Figure  7.10.  Velocity  profile 
at  location  x/H=4.1 


192 


1  . 

RE-  . 1700E+5 

IMAX-  41 

.  73 

s 

X  N. 

JMAX-  14 

Y 

N 

LOC: 

.  25 

\  K 

I-  31 

y. 

X-  11.17 

— / 

ISTEP-  lOO 

qne-sc>le 

U  • 

0 

• 

1 

5 

l  . 

,  5 

HOR . VEL 


Figure  7.11.  Velocity  profile 
at  location  x/H=7.1 


193 


Figure  7.12.  Turbulent  kinetic 
energy  profile  at  x/H=4.1 
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Figure  7.13.  Turbulent  kinetic 
energy  profile  at  x/H=7.1 
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that  the  two-scale  turbulent  model  indeed  predicts  better  flow  field  for 
the  two  dimensional  separation  flow  than  the  one-scale  model.  More  work 
is  needed  to  verify  the  prediction  capability  of  the  two-scale 
turbulence  model  in  other  flow  configurations.  Grid  dependence  studies, 
as  well  as  improvement  in  the  wall  function,  are  necessary  for  better 
prediction  of  turbulent  quantities. 
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CHAPTER  VIII 

CONCLUSIONS  AND  RECOMMENDATIONS 

8.1  Conclusions 

In  this  investigation,  a  new  two-scale  turbulence  model  is  developed. 
The  two  turbulent  scales  are  based  on  a  large  energy  containing  scale 
(k,s)  and  a  small  energy  dissipating  scale  (v,e).  The  (k,e)  scale  has 
been  used  in  previous  turbulence  models.  The  (v,e)  scale,  which  is  known 
as  the  Kolmogorov  scale,  is  used  in  the  present  investigation  to  model 
the  destruction  term  of  the  z -equation. 

The  two-scale  turbulence  model  shows  that  the  e -equation  needs  to 
have  the  influence  of  viscosity  since  viscosity  is  the  main  cause  of 
dissipation.  In  general,  the  two-scale  turbulence  model  predicts  a 
lower  kinetic  energy  than  the  one-scale  model. 

Calculations  of  free  shear  flows  in  Chapters  V  and  VI  and 
recirculating  flows  in  Chapter  VII  indicate  that  the  two-scale  k-e 
turbulence  model  gives  significant  improvement  over  the  one-scale 
turbulence  model.  It  is  important  to  point  out  that  unlike  the  one- 
scale  turbulence  model,  the  two-scale  turbulence  model  does  not  require 
modification  of  turbulent  constants  in  predicting  plane  or  round 
jets, mixing  layer  flows,  plane  and  round  wakes,  buoyant  jets,  flow  past 
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a  backward  facing  step  and  flow  past  a  channel  obstacle.  Therefore,  the 
predictability  of  the  two-scale  turbulence  model  is  demonstrated. 

Since  the  present  study  relied  on  the  existing  solution  techniques, 
it  should  be  mentioned  that  proper  understanding  of  the  numerical  scheme 
is  important  in  solving  the  flow  problem.  During  the  study  it  was  found 
that  it  is  necessary  to  choose  a  suitable  marching  step  in  the  modified 
Patankar  -  Spalding  algorithm  in  order  to  obtain  an  accurately  converged 
solution.  Also,  it  should  be  remarked  that  the  initial  conditions  affect 
the  solution  to  some  extent  particularly  near  the  inlet  or  initial  zone. 
Hence,  it  was  necessary  to  carefully  evaluate  experimental  data  before 
comparisons  between  either  various  experiments  or  experiment  and 
prediction  are  made. 

Further,  during  the  course  of  this  study,  it  was  found  that  though 
experimental  data  was  abundant,  many  sets  of  data  were  not  complete  for 
most  flows.  For  example,  numerous  results  were  given  at  various  sections 
of  the  flow  field,  but  initial  conditions  particularly  for  turbulent 
quantities  were  not  mentioned.  In  some  instances,  only  one  component  of 
normal  stress  was  available  for  comparison  instead  of  the  turbulent 
kinetic  energy  and  shear  stress.  Finally,  since  the  main  difference 
between  the  two-scale  and  the  one-scale  turbulence  model  is  in  the  t - 
equation,  it  would  be  of  interest  to  compare  the  dissipation  of 
turbulent  kinetic  energy  predicted  by  the  two  models.  However,  there 
were  almost  no  experimental  measurements  of  the  rate  of  dissipation  of 
turbulent  kinetic  energy. 
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8.2  Recommendation  for  future  srudv 
From  the  present  investigation  several  issues  still  require  further 
study  or  clarification.  They  are  explained  below. 

1.  In  chapter  II,  it  was  mentioned  that  the  turbulent  structure  in  a 
flow  is  a  function  of  the  Reynolds  number,  which  is  normally 
based  on  the  characteristic  velocity  and  length  of  the  problem. 
However,  in  most  experimental  investigations  it  is  the  turbulent 
Reynolds  number  based  on  the  Taylor  microscale  that  is  used.  The 
correlation  between  the  problem  or  mean  Reynolds  number  and 
Taylor  microscale  Reynolds  numbers  is  not  known.  Further  study  is 
necessary  in  turbulent  spectral  analysis  to  obtain  such  a 
relation.  It  is  worthwhile  to  see  how  the  turbulent  structure 
changes  directly  with  a  change  in  the  Reynolds  number  of  the 
problem.  This  may  verify  the  validity  of  the  turbulent  model  for 
a  range  of  Reynolds  numbers. 

2.  In  chapter  II  it  was  mentioned  that  the  spectral  analysis  was 
done  for  isotropic  flows  only.  At  present,  there  is  no  work 
available  for  analysis  of  nonisotropic  flows.  Such  an 
investigation,  if  possible,  could  help  in  understanding  the 
energy  transfer  process  and  provide  a  better  turbulence  model. 

3.  Various  aspects  of  the  computation  for  turbulent  flows  need  to  be 
considered  in  order  to  improve  the  two-scale  turbulence  model. 

For  example,  most  of  the  results  shown  in  Chapters  V,  VI  and  VII 
are  either  for  mean  velocity  profile,  centerline  velocity, 
turbulent  kinetic  energy  and  shear  stress.  It  is  desirable  to 
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further  investigate  the  prediction  of  the  two-scale  model  for 
other  flow  parameters  such  as  entrainment,  normal  stresses  and 
for  flows  involving  secondary  motion  and  three  dimensional 
configuration. 

4.  Further,  the  model  needs  to  be  extended  to  compressible  flows  and 
other  flows  with  strong  curvature. 
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APPENDIX  A 

COMPUTER  PROGRAM  GENMIX 


INSERT  SYSCOM>ERRD . INS . FTN 
INSERT  SYSCOM>KEYS . INS . FTN 
INSERT  SYSCOM>KEYS . INS . FTN 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

REAL* 8  LAB 

DIMENSION  OUT (83) ,DIFK(83) ,DIST(83) ,CONK(83) ,F4V(83) 

1 . DXRME (40) , XUU(40) ,YYH(40) ,YYHT(40) .XLONG(IOO) ,FLONG(8 , 100) 

2  ,FEI ( 100) 

DIMENSION  XPLOT(IOO) , YPLOT( 10 , 100) ,YAXES( 20) .SYMBOL (20) 

DIMENSION  FLUX ( 7 ) ,DFI(7) ,DFE(7) ,AJID(7) ,AJED(7) 

DIMENSION  YDT(83) ,UDT(83) ,TDT(83) ,CDT(83) ,FKDT(83) ,FEMDT(83) 
DIMENSION  FBI(7) ,FBE(7) ,EPS(83) ,FK3E2(83) ,FK4E3(83) 

DIMENSION  SYMBL1 (5 ) , SYMBL2 (4) 

COMMON/COMA/ A (83) ,AJE(7) ,AJI(7) ,B(83) ,C(83) ,CSALFA,D(83) ,DPDX(83) , 

1  DX,EMU(83) ,F(7 , 83) ,FS (5 , 83) , IAX, IEND , IFIN, INDE (7) , INDI (7) , IOUT, 

2  I STEP , ITEST , IUTRAP , JS , JSW , JV , JY , KEX , KIN , KRAD , N , ND2 , NF , NOVEL , NP 1 , 

3  NP2 ,NP3 ,OM(83) ,OMD(83) ,P(83) , PEI , PR(7) , PREF(7 , 83) , PSIE , PSII ,R(83) 
4 ,RHO(83) ,RME ,RMI ,RU(83) , SD(7 , 83) , SU(7 , 83) ,TAUE ,TAUI ,U(83) ,XD,XU, 
5Y(83) ,YE,YI ,ENU(83) ,NDEQ , BPI ,BPE ,DK1 (83) ,DK2(83) ,EDK1(83) ,EDK2(83) 
6 .US (83) .FACTOR 

COMMON/ COMB/ ARRCON , EWALL , H , HFU , INERT , MASSTR , MODEL , OXDFU , PREEXP , 

1  PRESS, UBAR.AK, RE, FR.ALMG.UFAC 

COMMON/ COMC/ENUT ( 83 ) ,ENUTDN(83) ,DUDY(83) ,DUDDY(83) ,DTDY(83) , 

1  DTDDY(83) ,PROK(83) ,BUPROK(83) ,ENUPR(83) ,PREFI(7) 

2, CDFN(83) ,PKDEP(83) ,CVFN(83) ,UV(83) , W(83) ,UT(83) , 

3  VT(83) ,TT(83) ,PTDET(83) ,PROT(83) ,DIFTT(83) , SUUK(83) , SUUE (83) 

4  ,SUUT(83) ,FUUK(83) ,FUUE(83) ,FUUT(83) 

COMMON/CONST/ IZT , CMUF , GDM , CDIS , Cl , C2 ,NCM , AL1 , ALD , CR1 , CRD , BUOY , 

1  C1T,C1K,C2K,NIB .ENULIM.NBUPRO .CT.CEl ,CE2 ,CT1 ,CC1 ,CC2 ,C2T, CSP ,CE 

2  , NPKDE , NPTDE , NALG , CAXIAL , PCLINR , LINEAR . LESSON , CE3 , CM2 , C2TM 
COMMON/ AUXL/RENO , VI SMIX , RHOA , TA , COEFEP , COEFED , EPSPK , EPSDT 

DATA  VI/  'TEST  l7,V2/'U7,V3/’F(l,I)7,V4/'F(2,I)7,V5/,F(3,I)7 

^DATA  V6 , V7 , V8 , V9 ,V10/ 'TEST2 ' , 'FS(1 , I) ' , 'FS(2,I)' , ’FS(3 , I) ’ , 'RHO(I) 

DATA  Vll , V12 , V13 , V14 , V15/ ' TEST3 ' , 'Y(I) ' , 'R(I) ' , 'RU(I) ' , ' TEST4 ' / 
DATA  V16  ,  V17  ,  V18  ,V19  ,  V20/  'EMU(I )  '  ,  '  OMEGA'  ,  'TEST5  '  ,  'R1  YS','VEL7 
DATA  V21 , V22 , V23  , V24 ,  V25  / 'TEMP' , 'T' , ’RHO' , 'KENGY' ,  77 
DATA  V26 ,V27 ,V28 ,V29 , V30/ 'DISSK' , 'D' , ’SU(I)' , ' ENUTDN ' , 'ENU'/ 

DATA  V31 , V32 , V33 , V34 , V35 / ' ENUT* , 'UV' , 'VV' , 'DIFTT' , 'DIST'/ 

DATA  V36,V37,V38,V39,V40/'UT' , 'VT' , 'DKDEP' , 'PTDET' , 'CD  FN'/ 

DATA  V4l.V42.V43 .V44.V45/ ' CV  FN' , 'TT-AGL' , 'TT-DEQ' , 'TT' , 'E '/ 
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DATA  V46 , V47 ,  V48 , V49 ,  V50/ ' DUDDY 1 , 1 DTDDY 1 , 'CONK* , 'DIFK' , 'PROK' / 

DATA  V51 , V52 , V53 , V54, V55/ 1 BUPROK1 , ' PREFTT 1 , ' PREFEP ' 

1, 'PREF-K' , ’PREF-T'/ 

DATA  SYMBL1/1HU, 1HT, 1HK, 1HD, 1HE/ 

DATA  SYMBL2/ 1HU, 1HT, 1HA, 1HF/ 

C 

CALL  SR(READ , ' INP1 ' , 4 , 7 ,TYPE , CODE) 

CALL  SR(WRIT, 'PSPl' ,4, 2, TYPE, CODE) 

CALL  SR (WRIT, 'PSP3' ,4, 3, TYPE, CODE) 

C 

C  . . - . . . . 

CHAPTER 11111111111111111  PARAMETERS  AND  CONTROL  INDICES  11111111111111 

c  .  NSTAT=  NO.  OF  STEPS  BETWEEN  OUTPUT  OF  SINGLE  VARIABLES. 

.  NPROF=  NO.  OF  STEPS  BETWEEN  OUTPUT  OF  ARRAY  VARIABLES. 

c  .  NPLOT=  NO.  OF  STEPS  BETWEEN  OUTPUT  OF  PLOT. 

c  .  IN  THIS  EXAMPLE,  PLOT  IS  CALLED  AT  END  OF  INTEGRATION  ONLY 

NSTAT=100 

NPR0F=1000 

NPL0T=20000 

C  .  MODEL=l  ONE -SCALE  K-E  MODEL 

C  .  MODEL=2  TWO -SCALE  K-E  MODEL 

MODEL=l 

c  .  IDAT=0=NO  DATA  READ  IN,  =1=  BUILT  IN  PROFILE 

IDAT=0 

C .  I LONG  =NO  OF  STEPS  BETWEEN  OUTPUT  IN  LONGITUDINAL  DIRECTION 

ILONG=300 

c  .  ISTTRAT=1=STRATIFCATION,=0=NO  STRATIFICATION 

I STRAT=0 

C  RATTD=RATE  OF  TD  VARIATION  IN  X  DIRECTION  (DEG. C/METER) -- 

RATTD=2 .0 

c .  LESSON  =  1  =  CE1*T0TAL  PROK  IN  EP  EQ. . — 

C  LESSON  =  2  =  CE1*SHEAR  PROK  AND  CE3  *  BUOY  PROD. IN  EP  EQ  - 

LESS0N=2 

C - NBY=1  WHEN  BUOYANT  FORCE  IS  INCLUDED - 

2  NBY=1 
7  KRAD=0 

C  . - . - . 

CHAPTER222222222222222222222222  GRID  AND  GEOMETRY  22222222222222222222 

—  —  —  —  —  «  —  —  —  —  —  —  -  —  —  —  —  — 

C  S.  I.  UNITS 

L=1 

C----N=NO  OF  GRID  POINTS 
N=40 
YN0Z=1 . 

XULAST=10 . OE+O 
LASTEP=20000 

C  IF(KRAD.EQ. 1 )LASTEP=3500 

C  IF(KRAD.EQ. 0)LASTEP=2100 

XU=1 . E-30 
XP=.99*XU 
X0UT=0 . 


o  o  o  on 
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XEND=0 . 

C .  FRA  =  PERCENTAGE  OF  FORWARD  MARCHING,  X  =  FRA  *  Y 

FRA=0 . 005 

C .  FACTOR  FOR  LINEARIZATION  OF  SOURCE  TERM  IN  E -EQUATION 

FACTOR3 1 . 0 
ULIM=. 01 
ENULIM=l.E+30 
TAN=0 . 0 1 
PEILIM=0 .01 
KIN=3 
KEX=2 
CSALFA=1 . 

-  R( 1)  ADJUSTMENT  MADE  JUST  BEFORE  CALL  STRIDE (1). 

. STRIDE4 . STRIDE4 . STRIDE4 . STRIDE4 

CALL  STRIDE (4) 

IAX=0 

IF (XEND . LE . 0 . )  IEND=0 
IF  (XOUT.LE.O.)  I0UT=0 

.  CHANGE  IEND ,  IAX  AND  IOUT,  IF  NECESSARY. 


CHAPTER333333333333333333333333  DEPENDENT  VARIABLES  SELECTION  33333333 


C  U(I)=VELOCITY 
C  F ( 1 , I )=STAGNATION  ENTHALPY 
C  F(2 , I)=CONCENTRATION  OF  FUEL 

C  F(3,I)=CONCENTRATION  OF  OXIDANT-OXDFU*F(2, I)=PHI 
C  FS(l,I)=CONCENTRATION  OF  OXIDANT 
C  FS ( 2 , I ) =TEMPERATURE 
C  FS ( 3 , I )=CONCENTRATION  OF  PRODUCT 
C  F(4,I)=KINETIC  ENERGY  OF  TURBULENCE 
C  F(5 , I)=DISSIPATION  RATE  OF  KINETIC  ENERGY 
C  F(6.I)=  THERMAL  ENERGY  OF  TURBULENCE  TT 
C  F(7 . I)=  DISS .RATE  OF  THERMAL  ENERGY 
C-  — -  FLONG(l ,L)=MAXIMUM  VELOCITY 
C— -  FLONG(2,L)=AXIAL  TEMPERATURE 

C -  FLONG(3,L)=ENTRAINMENT  COEFFICIENT 

C -  FEI (L)=BOUYANT  FLUX  INTEGRAL 

C  . NDEQ=0=ENTHALPY  AND  CONCENTRATION  EQUATION  NOT  SOLVED 

C  . NDEQ=1=ENTHALPY  EQUATION  SOLVED 

C  . NDEQ=2=ENTHALPY  AND  CONCENTRATION  EQUATION  SOLVED 

NDEQ=1 

NDEQP1=NDEQ+1 

NEQ=3+NDEQ 

C  ----NF=NO  OF  DEPENDENT  VARIABLES  TO  BE  SOLVED 
NF=6 

C - NDX=NO  OF  X  INTERVAL  USED  IN  CALCULATE  ENTRAINMENT  COEFF. 

NDX=40 

NDXM1=NDX-1 

C  . . . 


CHAPTER44444444444444444444444444444444444444C0NSTANTS  4444444444444444 
C  CHAPTER  4A . MATERIAL  CONSTANTS 


o  n 
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C  S . I .  UNITS 

C . GAS  CONSTAT  IN  JOULE/KILOMOLE/DEG.  K  . . 

GASCON=8300 . 

C - SPECIFIC  HEAT  AT  CONST.  PRESS.  IN  JOULE/KG/DEG . K - 

CMIX=1100. 

C . MOLECULAR  WEIGHT  IN  KG/KILOMOLE . 


WMIX=29. 

WA=29 . 

WD=29. 

GAMMA=CMIX/ (CMIX-GASCON/WMIX) 
V>rDGSCN=WMIX/GASCON 
VISMIX=1 . OE-7 


PREEXP=1 . 

C 

IZT=20 

C 

C  . CONSTANTS  FOR  TURBULENCE  MODEL 

CV=. 475 

C . CONSTANTS  FOR  UT  VT  EQS  - . 

C1T=3 . 2 
C2T=.5. 

C2TM=C2T 

C . CONSTANTS  FOR  UV  UU  W  EQS . 


CC1=2 . 8 
CC2=0 . 47 
CM2=CC2 
CSP=0 . 9 

IF (MODEL. EQ.l)  CSP=0.225 
CST=1 . 6 


C . CONSTANT  FOR  TT  EQUATION . ' 

CT1=0 . 13 
CT=I . 25 

C . CONSTANTS  FOR  TURBULENT  KINEMATIC  VISCOSITY 


CD=( 1 . -CC2)*CV/CC1 
PREFI (4)=1 . 

PREFI (5 )=1 . 3 
PREFI (6)=CD/CT1 
PREFI (7)=1 . 

CDIS=CD 

42  DO  40  J=1 , 3 
PR( J)=. 7 

40  PREFI ( J)=( 1 . -CC2)*C1T/CC1 
BUOY=FLOAT(NBY)*9 .81 

—  CONSTANTS  FOR  EP  EQ  -- 
CE=2 . 00 
CE1=17 . 50 
CE2=18 . 90 
CE3=CE1 

IF (MODEL. EQ. 1)  CE=0 . 15 
IF (MODEL. EQ. 1)  CE1=1.435 


o  n 
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IF (MODEL. EQ. 1)  CE2=1.92 
C0EFEP=CE1 
C0EFED=CE2 
WRITE (6, 45) 

WRITE (6, 43)  CSP, CC1, 002,011,021, CM2 ,  CE1  ,CE2 ,CE3 ,CE , CT1 , CT, 

1  FRA.C2TM 

45  FORMAT (20X,’ .  CONSTANTS  FOR  TURBULENCE  MODEL 

43  F0RMAT(4X, 'CSP  =' ,F10 . 5 ,5X, ' CC1  =' ,F10.5 ,5X, 'CC2  =  \F10.5,/ 

1  4X , ' C IT  =’ .F10.5.5X, 'C2T  = ' ,F10 . 5 ,5X, ' CM2  =’,F10.5,/ 

2  4X,  '  CE1  ='  ,F10 . 5 ,5X,  'CE2  = '  ,F10 . 5 , 5X,  ’  CE3  =',F10.5,/ 

3  4X , ' CE  =' ,F10 . 5 ,5X, ' CT1  =' ,F10 .5 ,5X, 'CT  =',F10.5,/ 

4  4X , ' FRA  =' ,F10.5,5X,’C2TM=' ,F10.5,5X/) 

C  . r  — 

CHAPTER55555 555555 555 555555 5555 5555555 55  INITIAL  CONDITIONS  5555555555 

C -  SPECIFY  RADIUS  , VELOCITY  ,TEMP .  ETC 

RA=0. 

RB=0 . 

RC=0 . 

RD=0 . 065 

DIAD=2.*RD 

UA=1 . 5 

UB=UA 

UC=UA 

UD=0 . 0 

TA=315. 

TD=300 . 5 

TDD=TD 

TB=TA 

TC=TA 

NEMU=1 

EMUI=1 . 

UREF=1 . 

YIN=RB 

TWALL=299. 

PRESS=1.E5 
PDGSCN=PRESS/ GASCON 
RHOA=PDGS  CN*WA / TA 
RHOD=PDGSCN*WD/TD 

C .  INITIAL  LONGITUDINAL  CONDITIONS  . 

DD=2.*RD 
XL0NG(1)=XU/DD 
FLONGCl, 1)=1. 

FL0NG(2,1)=1. 

FL0NG(3 , 1)=0 . 5 
FEI ( 1)=0 . 

--  CALCULATE  RE.FR  GR  RI  NUMBERS  . 

SBG=9.81*RATTD/TD 
BGDGD=RD*SBG/(9 . 81*(TA-TD)/TD) 

GUDGD=SBG*UA**2 . / ( (9 . 81*(TA-TD) /TD)**2 . ) 

REN0=(RH0A*UA*2*RD)/ (VISMIX*DSQRT(TA) ) 


non  ci  o 
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FRN0S=(UA**2*TA) / (9 .81*2*RD*(TA-TD)) 

FRN0=D  SQRT ( FRNOS ) 

GRN0=(9 .81*(TA-TD)*(2*RD)**3)/(TA*(VISMIX*DSQRT(TA) /RHOA) •■••2) 
RICHN=DSQRT (GRNO ) /RENO 
C 

WRITE (6,48)  RENO , FRNOS , FRNO , GRNO , RICHN , BGDGD , GUDGD , SBG 
48  FORMAT (/4X, 'RENO  =' ,E13 . 6 , 2X, ' FRN0S=' ,E13 . 6 , 2X, ’ FRNO  =',E13.6,/ 

1  4X , 1 GRNO  =' ,E13 . 6 ,2X, ' RICHN=' ,E13 . 6 , 2X, 1 BGDGD=' ,E13 . 6 , / 

2  4X, ' GUDGD=' ,E13 . 6 ,2X, ' SBG  =’,E13.6/) 

. CREATION  OF  INTIAL  PROFILES 

Y (NP3)=RD-RB 
EXPY=1 . 

UIN=UA 
UEN=UD 
TMPI=TA 
TMPE=TD 
C1I=1. 

C1E=0 . 

FATT=0 . 10 
FAET=0 . 10 
FA1=FATT 
FA2=FAET 
DO  52  J=4,NF 
FBI ( J)=0 . 

52  FBE ( J)=0 . 


INITIALIZE  SOME  OF  THE  PARAMETERS 

DO  53  1=1, NP3 
UT(I)=0 . 0 
VT(I)=0 .0 
UV(I)=0.0 
W(I)=0.0 
F4V(I)=0 . 0 
53  CONTINUE 
Y(1)=0. 

Y(2)=0 . 

YEPLS=0 . 0 
RME=0 . 0 
U1DUOL=0.0 
T1DT0L=0 . 0 
SB=1 . 0 
SA=1 . 0 
F(3,l)=1.0 
FE=1 . 0 

Y(NP2)=Y(NP3) 

U(1)=UIN 
U(2)=UIN 
U(NP3)=UEN 
U(NP2)=UEN+1 . OE-6 


non 
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UHB=U(1)-U(NP3) 

UHBS=UHB**2 
FS (2 , 1)=TMPI 
FS  (2 , 2)=TMPI 
FS(2,NP3)=TMPE 
FS (2 ,NP2)=TMPE 
F  (2, 1)=C1I 
F(2,2)=C1I 
F(2,NP3)=C1E 
F(2,NP2)=C1E 
DO  513  J=4,NF 

513  F ( J,NP3)=FBE ( J) 

IF (KIN. NE. 2)  GO  TO  510 
DO  514  J=4,NF 

514  F(J,1)=FBI(J) 

GO  TO  503 

510  F(4,1)=FA1*UHBS 

F (5 , 1)=0 . 10*F(4 , 1  )** l . 5 / ( 0 . 09*Y (NP3 ) ) 

F(6 , 1)=FATT*(FS (2 , 1) -FS (2 ,NP3) )**2 . 

F(7,l)=0. 10*F(6,1)**1.5/(0.09*Y(NP3)) 

503  DO  516  J=4,NF 
F(J,2)=F(J,1) 

516  F ( J,NP2)=F ( J,NP3) 

. PROFILES 

DYAF=1 . /FLOAT (N)**EXPY 
DO  520  1=2, NP2 

Y ( I ) = ( 1 . -DYAF*FLOAT (NP2 - I ) **EXPY) *Y (NP3 ) 

YP2=Y(I)**2. 

U(I)=U(NP3)+(U(1)-U(NP3))*DEXP(-2.8*YP2) 

FS (2 , I )=TMPE+(TMPI -TMPE )*DEXP ( -2 . 8*YP2/ 1.44) 
F(2,I)=C1E+(C1I-C1E)*DEXP(-2.8*YP2) 

F(3 , I)=l . -F(2 , I) 

F (4 , I )=UHBS*FA1*DEXP ( - 1 . 7*YP2) 

IF (KIN . EQ . 2 . AND . KEX.EQ . 2)  F(4, 1)=FA1*UHBS*(1. -( (Y(I)-Y(NP3)/2. ) 
1  /(Y(NP3)/2.))**2) 

F(5,I)«0.10*F(4,I)**1.5/(0.09*Y(NP3)) 
F(6,I)=F(6,1)*DEXP(-1.7*YP2) 

F(7,I)=0. 10*F(6,I)**1.5/(0.09*Y(NP3)) 

520  CONTINUE 
DO  521  1=1 ,NP3 
RHO(I)=PDGSCN/FS(2,I)*WMIX 

521  F(1,I)=CMIX*FS(2,I)+.5*U(I)*U(I)+F(4,I) 

IF (KIN. EQ. 3)  YNOZ=2 . *Y (NP3 ) 

IF(KIN.EQ.2)  YNOZ=2.*YIN 

C . CALCULATE  INITIAL  CD  FUNCTION  FOR  ENU(I)  . 

DO  525  1=1, NP3 
CDFN(I)=0 . 07 
525  CVFN(I)=0 .47 

WRITE (6 , 526)  PKDEP(3) ,CDFN(3) ,CVFN(3) ,N 
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526  FORMAT (/ ' PKDEP=' ,E13 . 6 , IX, ' CDFN  =' ,E13 . 6 , IX, ’ CVFN  =',E13.6,1X, 
l'N  =' , 15 ) 

DO  527  J=1 ,NF 
DO  527  1=1, NP3 

527  PREF(J, I)=l . 

C  CALCULATE  OMEGA  AND  STREAM  QUANTITIES  AND  WRITE 
501  YMDP=1 . 

DO  540  1=3, NP2 

IF(KRAD.EQ. 1)  YMDP=.5*(Y(I-1)+Y(I) )+YIN 
RUA= . 5* (RHO ( I ) *U ( I ) +RHO (I-1)*U(I-1)) 

540  0M(I)=0M(I-1)+YMDP*RUA*(Y(I) -Y(I-l) ) 

PEI=0M(NP2) 

C- --.--OMEGA , 0M(I ),  IS  MADE  DIMENSIONLESS  HERE  . - 

DO  541  1=3, NP2 

541  OM(I)=OM(I)/PEI 
PSII=YIN*U(1)*RH0(1) 

IFCKRAD.EQ. 1)  PSII=PSII*YIN/2 . 

PSIE=PSII+PEI 

C . SET  UP  INITIAL  SPREAD  PARAMETER  FOR  U  AND  T 

F1A=CMIX*TA+ . 5*UA**2 
F1D=CMIX*TD+ . 5*UD**2 
F2A=1 . 

F2D=0 . 

F3A=1 . -F2A 
F3D=1 . -F3A 
DYHA=0 . 

UBARDL=1 . 

DYHAV=1. 

YHA=0 . 

YHALS=0 . 25 
YELS=0 . 

DYHAT=0 . 

TBARDL=1. 

DYHAVT=1. 

YHAT=0 . 

YHATLS=0 . 25 
DO  549  K=1 ,NDX 
YYH(K)=0 . 

YYHT(K)=0. 

XUU(K)=0 . 

549  CONTINUE 

XUU(NDX)=XU 

C . SET  UP  INITIAL  ENTRAINMENT  VELOCITY . 

C  .  RTBDVB  IS  RATIO  OF  TEMP.  TO  VEL. LAYER  AT  BOUNDY,USE  IN 

c  .  ENTRAIN. CONTROL 

RTBDVB=1 .4 
ENTV=0 . 

ENTVLS=0 . 

DENTV=0 . 

YELS=RD 

IF(KIN.EQ.3)  U1V=U(1) 


o  n 
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FLOA=PSII 

IF(RC.LE.RB)  GO  TO  560 
DO  550  1=1, NP3 
YAB=Y ( I )+YIN 
IF(YAB.LT.RC)  GO  TO  550 
IC=I 

GO  TO  551 

550  CONTINUE 

551  FLOB=PEI*OM(IC) 

FLOC=PEI -FLOB 
ICM1=IC-1 
F2B=0 . 

DO  552  1=2 , ICM1 
OMDP=OM(I+l)-OM(I) 

TB=TB+(FS (2 , I)+FS (2 , 1+1) )*OMDP 

552  F2B=F2B+(F(2 , I)+F(2 , 1+1) )*OMDP 
TB=. 5*TB 

F2B=.5*F2B 
F3B=1 . -F2B 
F2C=0 . 

DO  553  I=IC,NP1 
OMDP=OM(I+l)-OM(I) 

TC=TC+(FS (2 , I)+FS (2 , 1+1) )*OMDP 

553  F2C=F2C+(F(2, I)+F (2 , 1+1) )*OMDP 
TC=. 5*TC 

F2C=.5*F2C 
F3C=1 . -F2C 
GO  TO  570 

560  FL0B=0. 

F2B=0 . 

F3B=0 . 

FLOC=PEI 
F2C=0 . 

DO  561  1=2 ,NP1 
OMDP=OM(I+l)-OM(I.) 

TC=TC+ (FS (2 , I)+FS (2 , 1+1) ) *OMDP 

561  F2C=F2C+(F(2,I)+F(2,I+l))*OMDP 
TC=.5*TC 

F2C=.5*F2C 
F3C=1 . -F2C 
570  FLOTOT=PEI 


CHAPTER6666666666666666666666666666666  THERMODYNAMIC  PROPERTIES  666666 
C  ir&iririririrfriririr!;  it  A"A' .V  A  it  START  OF  MAIN  LOOP 

60  CONTINUE 
C 

603  CONTINUE 


C 


IF(KIN.NE.2)  GO  TO  630 
IF(KRAD.EQ.l)  GO  TO  631 


ADJUST  R(l)  AND  YIN 


non  nnnnn  non 
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YIN=PSII/(U(1)*RH0(1)) 

GO  TO  630 

631  YIN=DSQRT (DABS (2 . *PSII/ (RH0(1)*U ( 1 ) ) ) ) 

R(1)=YIN 
630  CONTINUE 
C 

c  . STRIDE  1 . STRIDE  1 . STRIDE  1 . STRIDE1 

CALL  STRIDE (1) 

. RHO-'-U  R  AND  Y  ARE  CALCULATED 

. CALCULATION  OF  CHARACTERISTIC  FLOW  WIDTH  AND  UGL 

UHB=U( 1) -U (NP3) 

UHB2=1./(UHB**2) 

UHA=. 5*UHB 
UHC=DABS (UHA) 

UHA=UHA+U ( NP 3 ) 

UHE=DABS (UHB ) / CDEXP ( 1 . ) ) 

IHE=3 

-  YEP=  Y  AT  U=  UMAX/EXP  . 

-  YHA=  Y  AT  U=  UMAX/2  . 

-  YHAT=  Y  AT  T=  TMAX/2.  . 

DO  619  1=3, NP2 

IF(DABS(U(I)-U(NP3)) .GT.UHE)  GO  TO  619 

IHE=I-1 

GO  TO  618 

619  CONTINUE 

618  YEP=Y(IHE)+(UHE-U(IHE)  )*(Y(IHE+1)  -  Y(IHE) )  /  (U(IHE+1)  -U(IHE)  ) 

IHA=2 

DO  620  1=3, NP2 

IF (DABS(U(I) -U(NP3) ) . GT.UHC)  GO  TO  620 

IHA=I-1 

GO  TO  621 

620  CONTINUE 

. CALCULATE  Y  HALF  (YHA)  . 

621  YHA=Y(IHA)+(UHA-U(IHA))*(Y(IHA+1)-Y(IHA))/(U(IHA+1)-U(IHA)) 
YHR=YHA+YIN 

THB=FS(2, 1)-FS(2,NP3) 

THB2=1./(THB**2) 

THA=.5*THB  ^ 

THC=DABS (THA) 

THA=THA+FS ( 2 , NP3 ) 

IHT=1 

DO  626  1=3, NP2 

IF (DABS (FS(2,I)-FS(2,NP3)). GT.THC)GO  TO  626 
IHT=I - 1 

IF(IHT  .GT.  0)  GOTO  627 
626  CONTINUE 


ooo  oon  ooo 
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627  YHAT=Y ( IHT)+ (THA-FS ( 2 , IHT) )* (Y ( IHT+1 ) -Y(IHT) ) / 

1  (FS (2 , IHT+1) -FS (2 , IHT)+. IE -5) 

DYHAT=( YHAT-YYHT ( 1 ) ) / (XUU(NDX) -XUU ( 1 ) ) 

D YHA= ( YHA - YYH ( 1 ) ) / (XUU (NDX) -XUU ( 1 ) ) 

DYHATL= ( YHAT - YHATLS ) / ( XU - XP ) 

YHDYHT=YHA / YHAT 
DYEP=(YEP-YEPLS)/ (XU-XP) 

-CALCULTE  ENTRAINMENT  COEFFICINT . 

DYE=- (YELS2-YE)/ (XUU(NDX-2) -XP) 

ANG=ATAN (DYHA ) 

ANG=C0S (ANG) 

ENTV=RME*ANG/ (RHO (NP3)*R (NP3 ) ) /RTBDVB 
ALFA=ENTV/U(1) 

ALF  ALO=RME*ANG/ (RH0(NP3)*R(NP3)*U( 1 ) ) /RTBDVB 
ALFAAX=RME*ANG/ (RHO(NP3)*U(l)*(Y(IHE)**KRAD) ) /RTBDVB 
DENTV= ( ENTV - ENTVLS ) / (XU-XP) 

IF (KIN . NE . 2)  GO  TO  610 

. IF  FLOW  IS  OF  SHEAR  LAYER  . 


YHM=YHA 
UDI=. 9*UHB 
UDIC=DABS (UDI ) 

UDI=UDI+U(NP3) 

UDE=.1*UHB 

UDEC=DABS(UDE) 

UDE=UDE+U (NP3 ) 

YDI=0. 

IDI=2 

DO  622  1=3, NP1 

IF (DABS (U(I) -U(NP3) ) .GT.UDIC)GO  TO  622 
IDI=I-1 

IF(IDI  .GT.  0)  GO  TO  623 

22  CONTINUE 

23  YDI=Y(IDI)+(UDI-U(IDI) )*(Y(IDI+1) -Y(IDI) )/ (U(IDI+1) -U(IDI ) ) 
IDE=1 

DO  624  1=3, NP1 
IREAL=NP3-I+1 

IF (DABS (U ( IREAL) -U (NP3 ) ) . LT . UDEC )  GO  TO  624 
IDE=IREAL 

IF (IDE . GT. 0)  GO  TO  625 

24  CONTINUE 

625  YDE=Y(IDE)+(UDE-U(IDE) )*(Y(IDE+1) -Y(IDE) ) /  (U(IDE+1) -U(IDE) ) 
YHA=YDE-YDI 


CALCULATE  RATE  OF  SPREAD 


610  CONTINUE 

IF (U(NP3) . LE . 0 . )  GO  TO  640 


nnnnooooo 
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IF (DABS (UHB/U(NP3) ) . LT. .02)  DYHA=DYHA*U (NP3 ) /UHB 
640  DYHAL0= ( YHA - YHALS ) / DX 
YHALS=YHA 
YHATLS=YHAT 
YELS=YE 
YEPLS=YEP 
ENTVLS=ENTV 
YYH(NDX)=YHA 
YYHT ( NDX ) = YHAT 
XUU(NDX)=XP 
DO  650  K=1 ,NDXM1 
YYH(K)=YYH(K+1) 

YYHT(K)=YYHT(K+1) 

XUU(K)=XUU(K+1) 

650  CONTINUE 
C 

CHAPTER7 777777777777777777777777777777777777777  FORWARD  STEP  777777777 
DX=FRA*Y(NP3) 

IF ( I STEP . LT . 5  0 )  DX=DX/FLOAT ( 5 1 - 1 STEP ) 

C 

IF(ISTEP.GE.IEND)  GO  TO  70 
IF(DX.LT.XEND-XU)  GO  TO  70 
DX=XEND-XU 
IEND=ISTEP+1 

70  IF(ISTEP.GE.IOUT)  GO  TO  71 
IF(DX.LT.XOUT-XU)  GO  TO  71 
DX=XOUT-XU 
IOUT=ISTEP+l 

71  IF (DX . GT . XULAST-XU)  DX=XULAST-XU 
IF (DX . GT . 0 . )  GO  TO  73 

IFIN=1 
GO  TO  1011 
C 

73  XD=XU+DX 

-THIS  IS  JUST  INITIAL  GENERAL  SETUP 

FURTHER  ADJUSTMENTS  TO  DX  ARE  MADE  IN  CHAPTERS  8  AND  9. 


HAPTER8 88 8 88 88 8 8 88 8 88 8 8 888 8888  ADJUST  LONGITUDINAL  CONDITIONS  8888888 

CHAPTER8A . BOUNDARY  CONDITIONS 

.  I  BOUNDARY 

.  SYMMETRY  AXIS 

80  KIN=3 
RMI=0 . 

YIN=0 . 

R(1)=0. 

PSII=0. 

c  . - . FREE 


C 


85  KEX=2 


UBAR 


o  o  o  n  o  o  no  no  on 
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89  UBAR=0 . 

TBAR=0 . 

DO  820  1=2, NP 1 

TB AR=TB AR+ ( FS ( 2 , I ) +FS ( 2 , 1+ 1 ) ) *0MD ( I ) 

820  UBAR=UBAR+ (U ( I ) +U ( 1+1 ) )*0MD ( I ) 

TBAR= . 5*TBAR 
UBAR=.5*UBAR 

IF (KIN. EQ. 2)  UBAR= (UB AR -UA ) *PE I / PS IE+UA 

IF (KIN.EQ. 2)TBAR=(TBAR-FS (2 ,NP3) )*PEI/PSIE+FS (2 ,NP3) 


CHAPTER999999999999999999999  TRANSPORT  AND  ENTRAINMENT  PROPERTIES 

.  LAMINAR  VISCOSITY  ACCORDING  TO  SQUARE -ROOT  FORMULA 

WITH  WEIGHTING  ACCORDING  TO  MASS  FRACTION. 

DO  98  1=1, NP3 
FSAB=DABS (FS (2 , I) ) 

98  EMU ( I ) =VI SMIX*DSQRT ( FS AB ) 

CONTINUE 
EMU ( 2 ) =EMU ( 1 ) 

EMU ( NP  2 ) =EMU ( NP 3 ) 


9999 


. AUX . AUX . AUX . AUX . . AUX 

CALL  AUX 

-  SOURCE  TERM  FOR  K  AND  EP  ALSO  UV, W, VT,UT,TT, ARE  CALCULATED 

.  ENTRAINMENT  CONTROL. 

IF (KIN. NE. 2)  GO  TO  94 

RAT=DABS((U(3)-U(l))/(U(NP3)-U(l)+l.E-30)) 

IF (RAT . LT . ULIM)  EM  U(2)=EM  U(2)*RAT/ULIM 
RMI=2.*EM  U(2) 

94  CONTINUE 

IFCKEX.NE.2)  GO  TO  97 

RAT=DABS ( (U (NP1) -U (NP3) )/ (U(NP3) -U( 1)+1 .E-30) ) 

RME=-2.*EM  U(NP1)*RTBDVB 

IF (RAT . LT . ULIM)  RME=RME* (RAT/ULIM)**2 

IF (RAT. LT. ULIM*. 5)  RME=0 . 

97  IF (XD . EQ . XEND . OR . XD . EQ . XOUT . OR . XD . EQ . XULAST . OR . IAX .  EQ .  I STEP+1 ) 

1  GO  TO  96 

.  LIMIT  ON  INCREMENT  IN  PEI. 

IF ( (DABS (RMI )+DABS (RME) )*DX. LT. PEI*PEILIM)  GO  TO  96 
DX=PEI*PEILIM/ (DABS (RMI )+DABS (RME ) ) 

XD=XU+DX 
96  CONTINUE 


. STRIDE2 . STRIDE2 . STRIDE2 . STRIDE2 

95  CALL  STRIDE (2) 

-  SET  UP  A  B  C  D  ARRAY  . 


CHAPTER  10  10  10  10  10  10  10  10  10  10  10  10  10  10  10  OUTPUT  10  10  10 

C . SET  UP  OUTPUT  FORMAT . 

1000  AN STAT=NSTAT 


o  n  o 
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ANPROF=NPROF 


ANTPLOT=NPLOT 


IF (I STEP . GT. 0)  GO  TO  106 

CHAPTER  10A  .  HEADINGS 

REY=2 . *R (NP3) *RHO ( 1 ) *UBAR/ EMU ( 1 ) 

EQRAT=0 . 0 


WRITE (6 , 1013)  KRAD, LESSON 
C 


1013  FORMAT (/  'KRAD  =’ ,13, 10X, ’LESSON*’ ,13, 10X/) 

WRITE (6, 1010)  UA , UB , UC , UD , TA , TB , TC , TD , RA , RB , RC , RD , 
1  XULAST , PRESS , PREEXP 
C 


1010  FORMAT (/ 4X, 'UA 

1  ,4X,'UD 

2  ,4X,'TC 

3  , 4X , ' RB 

4  , 4X , ' XULAST= 
LAB=V17 


,E13 . 6 ,5X, 'UB 
,E13 . 6 ,5X, 'TA 
,E13 . 6 ,5X, 'TD 
.E13.6.5X, 'RC 
,E13 . 6 , 5X, 'PRESS 


,E13.6,5X, 'UC 
,E13.6,5X, 'TB 
,E13 . 6 ,5X, ' RA 
,E13 . 6 , 5X, ' RD 
,E13.6,5X, 'PREEXP* 


WRITE (6, 100)  LAB, (OM(I) ,I=1,NP3) 
PRESS 1=PRESS 
106  CONTINUE 


,E13 . 6/ 
,E13 . 6/ 
,E13 . 6/ 
,E13 . 6/ 
,E13 . 6/ ) 


IPRINT=0  GIVES  NO  OUTPUT,  =1  GIVES  SINGLE  VARIABLES  ONLY, 

=2  GIVES  BOTH  SINGLE  AND  ARRAY  (PROFILE)  VARIABLES. 

1011  IPRINT=0 

IF(FLOAT(ISTEP/NSTAT) .EQ. FLOAT (I STEP) /ANSTAT)  IPRINT=1 
IF (FLOAT (ISTEP/NPROF) . EQ. FLOAT (ISTEP)/ANPROF)  IPRINT=2 
IF(ISTEP.EQ. IEND.OR. ISTEP.EQ. IAX.OR. ISTEP.EQ. IOUT)  IPRINT=2 
IF(ITEST.NE.O.OR. IFIN.NE.O)  IPRINT=2 

c  .  THE  NEXT  STATEMENT  WOULD  BE  USED  FOR  A  TYPICAL  PLOT  CONTROL 

IF (FLOAT(ISTEP/NPLOT) . EQ. FLOAT (ISTEP)/ANP LOT)  IPRINT=3 

c  .  THE  NEXT  STATEMENT  PROVIDES  A  PLOT  JUST  PRIOR  TO  TERMINATION 

IF (XU . GE . XULAST . OR . IFIN . NE . 0 . OR . ISTEP . EQ . LASTEP)  IPRINT=3 
C 

1015  ALF AX=RME * ANG/ (RHO(IHE)*UBAR*(YEPLS**KRAD) ) 

C 

CHAPTER  IOC  . SINGLE  STATION  VARIABLES. 

1200  IF(IPRINT.EQ.O)  GO  TO  110 
UBAR=0 . 

DO  1021  1=2, NP1 

1021  UBAR=UBAR+0MD(I)*(U(I)+U(I+1) ) 

UBAR= . 5*UBAR 
UBARLS=UBARDL 
UBARDL=(UBAR-U(NP3) )/UHB 
DUBAR=DABS (UBARDL-UBARLS ) /UBARDL 
DDYHA=DABS ( (DYHA-DYHAV) ) 

DYHAV=DYHA 

U 1DU0=U ( 1 ) / (UA-U (NP3 ) ) 

T 1DT0= ( FS ( 2 , 1 ) -TDD ) / (TA -TDD ) 

—  TEST  FOR  DEVELOPED  FLOW  . 


C 
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IFCDUBAR . LT . 1 . E -3 . AND . DDYHA . LT . 1 . E-4 . AND . ISTEP . GT . 4000 )  IFIN=2 
IFCIFIN.EQ.2)  IPRINT=3 

WRITE (6,1030)  ISTEP , KRAD , KIN , KEX , FRNOS , PS I I , PS IE , U1DU0 , 

1  RMI , RME , PE I , U ( 1 ) , YHA , DYHA , UBARDL , T1DTO 

1030  FORMAT(////' ISTEP  ='  19  ,5X,  'KRAD  =',I9,5X, 

1  'KIN  =  ' ,I9,5X, 'KEX  =',I9,5X,/ 

1  'FRNOS  =' ,E9.3,5X, 'PSII  =' ,E9 . 3 ,5X, ' PSIE  =',E9.3,/ 

2  'U1DUO  =' ,E9 . 3 ,5X, 'RMI  =' ,E9 . 3 ,5X, 'RME  =’,E9.3,/ 

3  'PEI  =' ,E9.3,5X, 'U(l)  =' ,E9 . 3 ,5X, ' YHA  =',E9.3,/ 

4  'DYHA  =' ,E9 . 3 ,5X, ' UBARDL=' ,E9 . 3 , 5X, 'T1DTO  =' ,E9 . 3 , /) 
XDD=XU/DIAD 

DU1DU0=CU1DU0L-U1DU0)/ (XU-XP) /DIAD 
DT1DT0= (T1DT0L-T1DT0 )/(XU-XP)/DIAD 
WRITE (6,1014)  XDD , DU1DU0 , DT1DT0 , DYEP 
1014  FORMAT ('XDD  = ' E9 . 3 , 5X , ' DU1DU0= ' , E9 . 3 , 5X , ' DT1DT0= ' , E9 . 3 , 

1  5X/DYEP  ='  ,E9 . 3  ,  / ) 

U1DU0L=U1DU0 

T1DT0L=T1DT0 

WRITE (6, 1036)  DUB AR, DDYHA, B PE ,FS (2 , 1) ,FS(2,NP3) , YHAT , DYHAT , RATTD 

1036  FORMAT('DUBAR  =' ,E9 . 3 ,5X , ' DDYHA  =' ,E9 . 3 ,5X, ' BPE  =',E9.3,5X, 

1  'FS2I  =' ,E9.3,5X,/'FS2E  =' ,E9 . 3 ,5X, ' YHAT  =',E9.3,5X, 

2  'DYHAT  =  ' ,E9.3,5X, 'RATTD  =' ,E9.3,/) 

WRITE  (6,1037)  ALFA, ENTV.DENTV, DYE, CE1,CE2 

1 , ALFALO , DYHALO , ALFAX , DYHATL , YHDYHT 

1037  FORMATC ALFA  =' ,E9 -3,5X, 'ENTV  =' ,E9 . 3,5X, 'DENTV  =',E9.3,5X, 

1  'DYE  =  ’ ,E9.3,5X,/*CE1  =' ,E9 . 3 , 5X, ' CE2  =',E9.3,5X, 

2  ' ALFALO= ' , E9 . 3 , 5X , ' DYHALO= ' , E9 . 3 , 5X , / ' ALFAX  =',E9.3,5X, 

3  'DYHATL=' ,E9 .3,5X, 'YHDYHT=' ,E9 .3,/) 

DO  1020  J=1,NF 

1020  FLUX(J)=0. 

DO  1035  1=2, NP1 
DO  1035  J=1 ,NF 

1035  FLUX( J)=FLUX( J)+OMD(I)*(F ( J, I)+F ( J, 1+1) ) 

UFLUX=PE I *UB AR 
DO  1022  J=1 ,NF 

1022  FLUX ( J) = . 5*PEI*FLUX ( J) 

C 

UREF=U (NP3)+1 . E-30 
•  RUREF=UREF*RHO(NP3) 

AEXD=0 . 

DO  1023  J=1,NF 

DFI(J)=F(J, 1) -F(J,NP3)+1 .E-30 

1023  DFE(J)=DFI(J)+F(J, 1)-F(J,NP3) 

UFLUX=UFLUX-PEI*U(NP3)+U ( 1)*PSII 
GO  TO  (1041,1042,1043),  NDEQP1 

1043  FLUX(3)=FLUX(3) -PSIE*F3D+F3A*PSII 
1042  FLUX ( 2 ) =FLUX ( 2 ) - PS IE*F2D+F2A*PS I I 
1041  FLUX(1)=FLUX(1)-PSIE*F1D+F1A*PSII 
PRESSD=PRESS/PRESS1- 1 . 

RELO=U ( 1 ) *2 . *YHA*RHO ( 1 ) / ( VI SMIX*DSQRT (FS ( 2 , 1 ) ) ) 
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C 

WRITE ( 6 , 1 0  3 1 ) XU , DX , UFLUX , RELO 

1031  FORMAT ('XU  = ' ,E9 . 3 , 5X, ’ DX  =' ,E9 . 3 ,5X, ’UFLUX  =’ ,E9 . 3 ,5X, 

1  'RELO  =' ,E9 . 3 , /) 

C 

TK25=.5*(F(4,2)+F(4 , 3) ) 

£25=. 5* (F (5,2)+F(5,3)) 

C 

1026  CONTINUE 
C 

CHAPTER  10D  . PROFILES  AND  OTHER  ARRAYS. 

IF(IPRINT.EQ.l)  GO  TO  110 
C 

CHAPTER  10DD  . - . KINETIC  ENERGY  BALANCE 

C 

DO  1087  1=3, NP1 

FK3E2 (I)=F (4 , I)**3/F (5 , I)**2 

FK4E3 (I)=F(4, I )**4/F(5 , I)**3 

DIV=U(1)**3/YHA 

PROK(I)=PROK(I)/RHO(I)/DIV 

BUPROK ( I ) =BUPROK ( I ) /RHO ( I ) /DIV 

DIFK(I)=(ENU(I)*(R(I )+R(I+l) )*(F (4, 1+1) -F(4, I) )/ 

1  ( (Y(I+1) -Y(I ) )*PREF (4,1)) -ENU(I-1)*(R(I- 1)+R(I ) )*(F(4,I)-F(4,I-1) 

2  )/((Y(I)-Y(I-l) )*PREF (4, I) ) )/ ( (Y(I+1) -Y(I-1))*R(I)) 
DIFK(I)=DIFK(I)/RHO(I)/DIV 

CONK(I)=- (RHO(I)*U(I)*(F (4, I) -F4V(I) ) / (XU-XP)+PEI*(SA+SB*OM(I) ) 

1  *(F (4,1+1) -F(4,I-l))/( (Y(I+1) -Y(I-1))*R(I))) 
CONK(I)=CONK(I)/RHO(I)/DIV 
DIVT=FS(2,1)**2.*U(1)/YHAT 
PROT ( I ) =PROT ( I ) /RHO ( I ) /D IVT 

DIFTT(I)=(ENU(I)*(R(I)+R(I+1))*(F(6,I+1)-F(6,I))/ 

1  ( (Y(I+1) -Y(I) )*PREF (6,1)) -ENU(I-1)*(R(I-1)+R(I) )*(F(6 , I) -F (6,1-1) 

2  )/ ( (Y(I) -Y(I-l) )*PREF(6 , I) ) ) / ( (Y(I+1) -Y(I-l) )*R(I) ) 

D IFTT  ( I ) =D IFTT ( I ) /RHO ( I ) /DIVT 
DIST(I)=CT*F(5,I)*F(6,I)/F(4,I)/DIVT 

1087  CONTINUE 

C .  10E  .  10E  .  10E  . 

CHAPTER  10E  OUTPUT  TRAVERSE  PROFILES 
1086  CONTINUE 
LAB=V19 
DIV=YHA 

DO  1095  1=1, NP3 
1095  OUT(I)=Y(I)/DIV 

WRITE (6, 100)  LAB ,R(1) , (OUT(I) , 1=2 ,NP3) ,Y(NP3) ,YHA 
XAXIS=V12 
DO  1085  1=1, NP3 
1085  XPLOT(I)=OUT(I) 

LAB=V2 
SUB=0 . 

DIV=1 . 

IF (KIN. EQ. 3)  SUB=U(NP3) 


IF(KIN.EQ.3)  DIV=U(1)-U(NP3)+1  E-30 
DO  1094  1=1, NP3 
1094  OUT(I)=(U(I)-SUB)/DIV 

WRITE (6, 100)  LAB ,U( 1) , (OUT(I) , 1=2 ,NP3) ,U(NP3),DIV 
NY=1 

YAXES (NY)=V20 
SYMBOL (NY)=SYMBL1 (NY) 

DO  1084  1=1, NP3 
1084  YPLOT (NY , I ) =OUT ( I ) 

IF(NDEQ.EQ.O)  GO  TO  1091 
C 

LAB=V21 
SUB=FS (2 ,NP3) 

DIV=1 . E-30+FS(2 , 1) -FS (2 ,NP3) 

DO  1093  1=1, NP3 
1093  OUT(I)=(FS (2 , I) -SUB) /DIV 

WRITE (6, 100)  LAB, FS (2,1), (OUT(I) ,I=2,NP3) ,FS (2 ,NP3) ,DIV 

NY=NY+1 

YAXES ( NY) =V21 

SYMBOL (NY)=SYMBL1 (NY) 

DO  1083  1=1, NP3 
1083  YPLOT (NY , I ) =OUT ( I ) 

C 

LAB=V23 

WRITE (6 , 100)  LAB , (RHQ ( I ) , 1=1 , NP3 ) 

C 

1091  CONTINUE 
9999  LAB=V24 

DIV=(U(1)-U(NP3))**2 
DO  1097  1=1, NP3 

1097  0UT(I)=F(4, I) /DIV 

WRITE(6 , 100)  LAB, (OUT(I) ,I=1,NP3) ,DIV 

NY=NY+1 

YAXES (NY)=V24 

SYMBOL (NY)=SYMBL1 (NY) 

DO  1098  1=1, NP3 

1098  YPLOT(NY, I)=OUT(I) 

C 

LAB=V26 

DIV=(U(1)-U(NP3))**3/YHA 
DO  1096  1=1, NP3 
1096  OUT(I)=F(5 , I) /DIV 

WRITE(6 , 100)  LAB , (OUT(I) , 1=1 ,NP3) ,DIV 

NY=NY+1 

YAXES (NY)=V26 

SYMBOL (NY) =SYMBL1 (NY) 

DO  1099  1=1, NP3 

1099  YPLOT (NY, I)=OUT(I) 

DITT=(FS (2 , 1) -FS (2 ,NP3) )**2 . 

IF(ISTEP.EQ.O)  GO  TO  1114 
LAB=V28 


WRITE (6, 100)  LAB , (US (I) , 1=1 ,NP3) 

LAB=V29 

WRITE (6, 100)  LAB , (ENUTDN(I) , 1=1, NP3) 
LAB=V30 

WRITE(6 , 100)  LAB, (ENU(I) ,1=1, NP3) 
LAB=V31 

WRITE (6 , 100)  LAB, (ENUT(I) ,1=1, NP3) 
LAB=V32 

DIV=(U(1)-U(NP3))**2 
DO  1071  1=1, NP3 
1071  OUT(I)=UV(I)/DIV 

WRITE(6 , 100)  LAB, (OUT(I) , 1=1, NP3) ,DIV 

IF (MODEL. NE. 4)  GO  TO  1120 

LAB=V33 

DO  1100  1=1, NP3 
1100  OUT(I)=W(I)/DIV 

WRITE  (6,100)  LAB  ,  (OUT(I)  ,  1=1  ,NTP3)  ,DIV 
LAB=V34 

WRITE (6 , 100) LAB , (DIFTT(I) , 1=1 ,NP3) 
LAB=V35 

WRITE (6 , 100) LAB , (DIST(I) , 1=1 ,NP3) ,DIVT 
DIUT=U( 1)*(FS(2 , 1) -FS (2 ,NP3) ) 

LAB=V36 

DO  1105  1=1, NP3 
1105  OUT(I)=UT(I)/DIUT 

WRITE ( 6 , 100 ) LAB , (OUT ( I ) , 1=1 , NP3 ) , DIUT 
LAB=V37 

DO  1110  1=1, NP3 
1110  OUT(I)=VT(I)/DIUT 

WRITE (6, 100)LAB, (OUT(I); 1=1, NP3), DIUT 
LAB=V38 

WRITE (6 , 100) LAB , (PKDEP(I) , 1=1 ,NP3) 
LAB=V39 

WRITE (6, 100) LAB, (PTDET(I) ,I=1,NP3) 
LAB=V40 

WRITE(6, 100)  LAB , (CDFN(I) , 1=1, NP3) 
LAB=V41 

WRITE(6 , 100)  LAB , (CVFN(I) , 1=1 ,NP3) 
LAB=V42 

DO  1115  1=1, NP3 
1115  OUT(I)=TT(I)/DITT 

WRITE ( 6 , 100 ) LAB , ( OUT ( I ) , 1=1 , NP3 ) , DITT 
1114  LAB=V43 

DO  1111  1=1, NP3 
1111  OUT(I)=F(6 , I) /DITT 

WRITE (6, 100) LAB,  (OUT(I) , 1=1, NP3) , DITT 
IF(ISTEP.EQ.O)  GO  TO  110 
NY=NY+1 
YAXES (NY)=V44 
SYMBOL (NY) =SYMBL1 (NY) 

DO  1113  1=1, NP3 
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1113  YPLOT ( NY , I ) =OUT ( I ) 
1120  CONTINUE 


C 


LAB=V46 

WRITE(6 , 100)  LAB , (DUDDY(I) , 1=1, NP3) 
LAB=V47 

WRITE (6 , 100)  LAB, (DTDDY(I) ,I=1,NP3) 
LAB=V48 

WRITE (6 , 100)  LAB , (CONK(I) , 1=1 ,NP3) 
LAB=V49 

WRITE(6 , 100)  LAB , (DIFK(I) , 1=1 ,NP3) 
LAB=V50 

WRITE (6 , 100)  LAB , (PROK(I) , 1=1 ,NP3) 
LAB=V51 

WRITE(6 , 100)  LAB, (BUPROK(I) , I=1,NP3) 
LAB=V5  2 

WRITE (6 , 100) LAB , (PREF (6 , I ) , 1=1 ,NP3 ) 
LAB=V53 

WRITE (6 , 100) LAB , (PREF (5 ,  I ) , 1=1 ,NP3) 
LAB=V54 

WRITE (6 , 100)  LAB, (PREF (4, I) ,1=1, NP3) 
LAB=V55 

WRITE(6, 100)  LAB , (PREF(1 , I) , 1=1, NP3) 


1009  CONTINUE 

IFCIPRINT.EQ.2)  GO  TO  110 
IF(ISTEP.EQ.O)  GO  TO  110 
WRITE (6, 1070)  XU , I STEP , KRAD , FRNOS , DYHALO , DYHATL 
1070  FORMATC'PLOT  AT  XU=' ,E9  3, IX, 'ISTEP=' ,15, IX, 'KRAD=' ,15, IX, 
1  ’FRNOS=' , E9 . 3 , IX , 1 DYHALO= ' ,E9.3,1X, ' DYHATL= 1 ,E9.3/) 

CALL  PLOTS (XPLOT , 43 , NP3 , XAXIS , YPLOT , 10 , NY , YAXES , SYMBOL) 
NPLOT2=2 

IF (NPLOT2 . NE . 1 )  GO  TO  110 
NY=1 


YAXES (NY)=V32 
SYMBOL (NY )=V64 
DO  1201  1=1, NP3 
1201  YPLOT (NY , I ) =UV ( I ) 
NY=NY+1 
YAXES (NY) =V33 
SYMBOL (NY) =V65 
DO  1210  1=1, NP3 
1210  YPLOT(NY,I)=W(I) 
NY=NY+1 
YAXES (NY)=V36 
SYMBOL (NY) =V66 
DO  1220  1=1, NP3 
1220  YPLOT (NY, I)=UT(I) 
NY=NY+1 
YAXES (NY)=V37 
LAB=V67 

DO  1230  1=1, NP3 


on  no 
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1230  YPLOT (NY , I ) =VT ( I ) 

CALL  PLOTS (XPLOT , 43 ,NP3 , XAXIS , YPLOT , 10 ,NY , YAXES , SYMBOL) 


CHAPTER  11  11  11  11  11  11  11  ll  ii  ii  ii  END  OF  MAIN  LOOP 

110  IF (I STEP . GE . LASTEP . OR . XU . GE . XULAST. OR. IFIN . NE . 0)  GO  TO  111 
XP=XU 

SA=RMI/PEI 
SB=(RME-RMI)/PEI 
DO  113  1=1, NP3 
113  F4V(I)=F(4, I) 

C 

LX=L*ILONG 

IF(ISTEP . NE . LX)GO  TO  114 
L=L+1 

XLONG ( L ) =XU / DD 
FLONG(l,L)=(U(l) -U(NP3) ) 

FLONG(2 ,L)=(FS (2,1) -FS(2 ,NP3) )/ (TA-TD) 

FLONG(3 ,L)=-ALFAX 
AFDTT=DABS (F(6 , 1)/DITT) 

FLONG ( 4 , L ) =DS  QRT ( AFDTT ) 

IF(NBY.EQ.O)  GO  TO  104 

FLONG ( 5 , L ) =FLONG ( 1 , L ) *XLONG ( L ) ** ( KRAD / 3 . ) 

FLONG ( 6 , L)=FL0NG ( 2 , L)*XLONG(L)**( ( 3 . +2 . *KRAD ) / 3 . ) 

GO  TO  105 
104  CONTINUE 

FLONG ( 5 , L ) =FLONG ( 1 , L ) *XLONG ( L ) ** ( ( 1 . +KRAD ) / 2 . ) 

FLONG ( 6 , L ) =FLONG ( 2 , L ) *XLONG ( L ) ** ( ( 1 . +KRAD ) / 2 . ) 

105  FEI (L)=FE 

FLONG (7 ,L)=F(4,1)/((U(1) -U(NP3) )**2) 

c  . STRIDE3 . STRIDE3 . STRIDE3 . STRIDE3 

114  CALL  STRIDE (3) 

IF(IFIN)  1011,60,111 

.  TERMINATION 

111  WRITE (6 , 112)  I STEP, LASTEP, XU, XULAST, IFIN, DDYHA 

112  FORMAT (23H  TERMINATED  WITH  ISTEP=,I5,8H  LASTEP=,I5, 

1  4H  XU=, 1PE11 . 3 , 8H  XULAST=,E11 . 3 , 6H  IFIN=,I3,7H  DDYHA= , E 1 1 . 3 ) 

.  HA  .  11A  .  11A  .  11A  . 

----  CHAPTER  11A  OUTPUT  LONGITUDINAL  PROFILES 
NLONG=l 

IF (NLONG . EQ . 0 ) GO  TO  120 
WRITE (6 , 103) 

DO  115  1=1, L 

115  WRITE(6, 102)1, XLONG(I) ,FL0NG(1 , I ) ,FLONG(2 , I) ,FLONG(3 , I) ,FEI (I) 

1  , FLONG (4, I) , FLONG (5 , I ) , FLONG (6,1), FLONG (7,1) 

102  F0RMAT(I4, 1P11E11.3) 

103  FORMAT ( 3X , ' L 1 ,4X, 'X/D' ,4X, 'U(1)/UA' ,2X, ' (TC-TD) / (TA-TD) ' , 

1  2X , ' -ALF AX' , 3X , ' FE  INTGL' , IX, 'DSQRT(TT)/ (TC-TA) ' , IX, 'UAXIAL*X' 

2  , IX , ' TAXIAL*X ' , IX , ' KENG ' ) 

XAXIS=V68 

IF (L . GT. 65)  L=65 


n  n 
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IF(L  .EQ.  65)  GO  TO  120 
DO  1310  1=1, L 
1310  XPLOT ( I ) =XLONG ( I ) 

NY=1 

DO  1330  J=1 , 3 
YAXES (NY)=V69 
DO  1320  1=1, L 

1320  YPLOT(NY, I)=FLONG( J, I) 

1330  NY=NY+1 

SYMBOL ( 1 )=SYMBL2 ( 1 ) 

SYMBOL(2)=SYMBL2 (2) 

SYMBOL ( 3 ) =S YMB L2 ( 3 ) 

SYMBOL (4 )=SYMBL2 (4) 

YAXES (NY)=V72 
DO  1340  1=1, L 
1340  YPLOT(NY, I )=FEI (I) 

1350  FORMAT( 1H1 , 35H  PLOT  LONGITUDINAL  VARIABLES  AT  Y=0) 

CALL  PLOTS (XPLOT , 65 , L,XAXIS , YPLOT , 10 , NY , YAXES , SYMBOL) 


CHAPTER  11B  CONTROL  FOR  CASES  CALCULATION 
120  CONTINUE 
150  DO  119  1=1, NP3 
YY=Y(I)/YHA 

UU=(U(I)-U(NP3))/(U(l)-U(NP3)+l.E-30) 
WRITE(7,8888)  YY,UU 
119  CONTINUE 

DO  129  1=1, NP3 
YY=Y(I)/YHA 

FF=F(4,I)/(U(1)-U(NP3))**2 
C  FF=F(4,I)/(U(1)**2) 

WRITE(7,8888)  YY,FF 
129  CONTINUE 

DO  139  1=1, NP3 
YY=Y(I)/YHA 

FUV=UV ( I ) / (U ( 1 ) -U (NP3 ) ) **2 
C  FUV=UV(I)/(U(1)**2) 

WRITE(7 ,8888)  YY.FUV 
139  CONTINUE 

DO  149  1=1, NP3 
YY=Y(I)/YHA 

FT=(FS (2 , I) -FS(2 ,NP3) )/ (FS (2 , 1) -FS (2 ,NP3) ) 
WRITE (7 , 8888)  YY,FT 
149  CONTINUE 

DO  159  1=1, L 

UUU=FLONG (1,1)/ (DSQRT (UA* (UA-UD ) ) ) 

WRITE (7 , 8888)  XLONG(I),UUU 
159  CONTINUE 
8888  FORMAT (  2E20.7) 

CALL  EXIT 

C  100  FORMAT ( 1H  , A8 , 1P11E10 . 3 , 8 (/9X, 11E10 . 3) ) 

100  FORMAT (//1H  , A8 , 1P5E14 . 7 , 8 (/9X.5E14 . 7) ) 
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C  100  FORMAT ( 1H  ,A8 , 1P11E10 . 3/ (9X, 11E10 . 3)/ (9X, 12E10 . 3) ) 

101  FORMAT ( 1H  ,A8, 11111) 

END 

SUBROUTINE  AUX 

C/ . SUBROUTINE  FOR  PROGRAM  GENMIX  4A 

IMPLICIT  DOUBLE  PRECISI0N(A-H,0-Z) 

REAL* 8  LAB 

DIMENSION  YMPIC83) ,UAV(83) ,Z(83) ,FLAV(83) ,RAVE(83) 

DIMENSION  DUDO(83) ,SC(83) ,SCV(83) ,GD(83) ,DS(2,43) , YEDGE (6) 
COMMON/COMA/A(83) ,AJE(7) ,AJI (7) ,B(83) ,C(83) ,CSALFA,D(83) ,DPDX(83) , 

1  DX,EMU(83) ,F(7 , 83) ,FS (5 , 83) , IAX, IEND , IFIN , INDE(7) , INDI (7) , I OUT, 

2  ISTEP , ITEST, IUTRAP , JS , JSW , JV , JY ,KEX , KIN ,KRAD ,N,ND2 ,NF .NOVEL, NP1 , 

3  NP2 ,NP3 ,OM(83) , OMD(83) ,P(83),PEI,PR(7), PREF (7 , 83) ,PSIE , PSII ,R(83) 
4,RH0(83) ,RME ,RMI ,RU(83),SD(7,83),SU(7,83) ,TAUE ,TAUI ,U(83) ,XD,XU, 
5Y(83)  ,YE  ,YI  ,ENU(83)  ,NDEQ,BPI  ,BPE,DK1(83)  ,DK2(83)  ,EDK1(83) ,EDK2(83) 
6, US (83) , FACTOR 

COMMON/ COMB / ARRCON , EVALL , H , HFO , INERT , MAS  STR , MODEL , OXDFU , PREEXP , 

1  PRESS, UBAR.AK, RE, FR,ALMG,UFAC 

C0MM0N/C0MC/ENUT(83) ,ENUTDN(83) ,DUDY(83) ,DUDDY(83) ,DTDY(83) , 

1  DTDDY(83) ,PROK(83) ,BUPROK(83) ,ENUPR(83) ,PREFI(7) 

2  ,CDFN(83) , PKDEP(83)  ,CVFN(83)  ,UV(83)  ,W(83)  ,UT(83)  , 

3  VT(83) ,TT(83) ,PTDET(83) ,PR0T(83) ,DIFTT(83) ,SUUK(83) ,SUUE(83) , 

1  SUUT(83) ,FUUK(83) ,FUUE(83) ,FUUT(83) 

COMMON/CONST/ IZT , CMUF , GDM , CDI S , Cl , C2 , NCM , AL1 , ALD , CR1 , CRD , BUOY , 

1  C IT , C IK , C2K , NIB , ENULIM , NBUPRO , CT , CE 1 , CE2 , CT1 , CC 1 , CC2 , C2T , CSP , CE 

2  , NPKDE , NPTDE , NALG , C AXIAL , PCLINR , LINEAR , LESSON , CE3 , CM2 , C2TM 
COMMON/ AUXL/RENO , VI SMIX , RHOA , TA , COEFEP , COEFED , EPSPK , EPSDT 

DATA  VI/ 'TEST  1 1 / , V2/ 'U'/ ,V3/'F(1,I)'/ , V4/ ,F(2,I),/,V5/'F(3,I)'/ 
DATA  V6,V7,V8,V9,V10/,TEST2’  ,  'FS(1,I)' ,  'FS  (2 , 1)  '  ,  '  FS  (3 , 1)  '  ,  'RHO(I) 

1  7 

DATA  Vll , V12 , V13 , V14 , V15/ 'TEST3 ' , 'Y(I) ’ , 'R(I) ’ , ’RU(I) ’ , ’TEST4 V 
DATA  V16,V17,V18,V19 ,V20/ 'EMU(I) ' , 'OMEGA' , 'TESTS  1 , 'R1  YS’.’VELV 
DATA  V21, V22 ,V23 ,V24,V25  / 'TEMP' , 'T' , 'RHO' , 'KENGY' , 'K'/ 

DATA  V26 ,  V27  ,  V28  ,V29  ,  V30/  'DISSK'  ,  '  D '  ,  'SU(I)  '  ,  '  ENUTDN'  ,  'ENUV 
DATA  V31,V32,V33,V34,V35/'ENUT'  ,  'UV'  ,  'W'  ,  'DIFTT'  ,  'DISTV 
DATA  V36 ,  V37  ,  V38  ,V39  ,  V40/  'UT'  ,  '  VT'  ,  'DKDEP '  ,  '  PTDET1  ,  'CD  FN7 
DATA  V41,V42,V43,V44,V45/ 'CV  FN' , ' TT-AGL' , 'TT-DEQ ' , ' TT' , ' E 7 
DATA  V46,V47,V48,V49,V50/ 'DUDDY' , ' DTDDY ' , 'CONK' , 'DIFK’ , 'PROK'/ 
DATA  V51 , V52 , V53 , V54 , V55/ ' BUPROK' , ' PREFTT ' , 'PREFEP' 

1, 'PREF-K' , ' PREF -T 7 
C 

C  SD(3 , I)  IS  USED  FOR  R(I )*(Y(I+1) -Y(I-l) ) 

C . 

C .  K  E  MODEL  . 

C . CALCULATE  V  AND  T  GRADIENT 

500  DO  550  1=2, NP2 
CDFN(I)=0 . 07 
550  CVFN(I)=0 .47 
600  DO  601  1=2, NP1 

YMPI (I)=Y(I+1)-Y(I) 

DUPI=U(I+1)-U(I) 


non 
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DUDY ( I ) =DUP I / YMP 1(1) 

DTDY(I)=(FS(2, 1+1) -FS(2, I) )/YMPI(I) 
UAV(I)=.5*(U(I)+U(I+1)) 

CDFN(I)=. 5*(CDFN(I+1)+CDFN( I ) ) 

601  DUDO ( I ) =DUP I / OMD ( I ) 

DUDY ( 1 ) =0 . 

DTDY ( 1 ) =0 . 

DUDY (NP2)=DUDY (NP1) 

DTDY (NP2)=DTDY (NP1) 

CVFN (NP3 )=CVFN (NP2 ) 

CDFN (NP3 )=CDFN (NP2 ) 

DO  630  I  =2 ,NP1 


C-- . CALCULATE  ,K=Z(I),  EPS ILON=FLAV (I) , RAVE=R  AVERAG - 

Z(I)=.5*(F(4,I)+F(4,I+1)) 

FLAV(I)=. 5*(F (5 , I)+F(5 , 1+1) ) 

RAVE (I)=. 5*(R(I)+R(I+1) ) 

C . CALCULATE  TURBULNT  VISCOSITY 


628  ENUT ( I )=0 . 5* (RHO ( I )+RHO ( 1+1 ) )*CDFN ( I )*Z ( I )**2/FLAV ( I ) 

ENUTDN ( I )=1 .  +  (2 .  /  (2 . *C1T- 1 .  +PKDAVE+CT* (PTDAVE - 1 . )  )*Z  ( I ) 

1/FLAV(I) )*(BUOY/FS (2 ,NP3) )*(DTDY(I)/DUDY(I ) )*( ( 1 . -CM2) / (1. -CC2)) 
IF(ENUTDN(I) .GT.ENULIM)  ENUTDN ( I )=ENULIM 
630  ENU(I)=ENUT(I)*ENUTDN(I) 

C . CALCULATE  REY.  STRESS  UV . 

ENU(1)=ENU(2) 

ENU(NP2)=ENU(NP1) 

ENU(NP3)=ENU(NP2) 

ENUT(1)=ENUT(2) 

ENUT(NP2)=ENUT(NP1) 

ENUT (NP3 )=ENUT (NP2 ) 

FLAV ( 1 ) =FLA V ( 2 ) 

FLAV (NP2 ) =FLAV (NP 1 ) 

Z(1)=Z(2) 

Z(NP2)=Z(NP1) 

ENUTDN(1)=1. 

ENUTDN ( NP 2 ) =ENUTDN ( NP 1 ) 

ENUTDN ( NP  3 ) =ENUTDN (NP2) 

UV(NP3)=0. 

UV(1)=0. 

UV (NP2)=(UV (NP1)+UV(NP3) )* . 5 


22222222222222222222222222  VISCOSITIES 

.  LAMINAR  VISCOSITIES  FOR  CELL  BOUNDARIES 

200  DO  23  1=2, NP1 
23  EMU(I)=.5*(EMU(I)+EMU(I+1) ) 


C  . - .  TURBULENT  CONTRIBUTION 

DO  20  1=2, NP1 

EMU ( I ) =EMU ( I ) +ENU ( I ) 

20  CONTINUE 

C  .  MODIFICATION  OF  EMU  ARRAY 

29  DO  24  1=2, NP1 

24  EM  U(I)=EMU(I)/ (Y(I+1) -Y(I) ) 


o  o  o  o  o  o 
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IF(KRAD.EQ.O)  GO  TO  25 
DO  26  1=2, NT  1 

26  EM  U(I)=EM  U(I)*.5*(R(I)+R(I+1)) 

25  CONTINUE 

C  . . - . - .  INITIAL  PREF  S. 

DO  230  1=1, NP2 

PREF (1,I)=(1.-CC2)*C1T/CC 1*ENUTDN ( I ) 

PREF(2,I)=1. 

PREF(3 , I)=l. 

PREF (4, I)=l . 0*( 1-CC2) / (CSP*CC1) 

PREF (5 ,I)=1 .  0*(1-CC2)/ fCE*CCl) 

PREF (6 , I)=2 . *ENU(I )*FLAV(I)/ (CTl*(RHO(I)+RHO( 1+1) )*Z(I)*Z(I)) 

230  PREF(7 , I)=l . 

DO  231  J=1,NF 

231  PREF(J,NP3)=PREF(J,NP1) 

DO  227  1=1, NP3 

227  ENUPR(I)=ENU(I ) /PREF (1,1) 

333333333333333333  SOURCES 
.  VELOCITY  U 

- CALCULATE  SOURCE  TERM  IN  M  EQ  . 

DO  308  1=3 ,NP1 

308  US(I)=S  D(3 , I)*(DPDX(I) -BUOY*RHO(I)*(FS (2 , I) -FS(2 ,NP3) ) /FS (2 ,NP3) ) 
US (2)=S  D(3,2)*(DPDX(2)-BU0Y*RH0(1)*(.25*(2.*FS(2,1)+ 

1  FS(2,2)+FS(2,3))/FS(2,NP3)-1.)) 

US (NP2)=S  D(3,NP2)*(DPDX(NP2)-BUOY*RHO(NP3)*(.25*(2.*FS(2,NP3)+ 

1  FS (2 ,NP2)+FS (2 ,NP1) )/FS (2 ,NP3) -1 . ) ) 

.  TO  CALCULATE  SOURCE  TERM  FOR  K  AND  EPS  EQS  - 

.  FIRST  COMPUTE  W,VT,TT,UT,UV,ETC-IN  CENTRAL  DIFFERENCE . 

. K  AND  EPS 

800  DO  801  1=2, NP2 

DTDDY ( I ) = . 5  * (DTDY ( I ) +DTDY ( I - 1 ) ) 

801  DUDDY(I)=.5*(DUDY(I)+DUDY(I-1)) 

DTDDY ( 1 ) =0 . 

DUDDY(1)=0 . 

F4KIN=. 5* ( . 5*(F(4,3)+F(4 , 2) )+F(4 , 1) ) 
F5KIN=.5*(.5*(F(5,3)+F(5,2))+F(5,1)) 

F(4,2)=F4KIN 

F6KIN=.5*( . 5*(F(6 , 3)+F(6 , 2) )+F(6 , 1) ) 

F (5 , 2)=F5KIN 
F(6 , 2)=F6KIN 

F (4 ,NP2)=0 .5*( . 5*(F(4 ,NP1)+F (4 ,NP2) )+F(4 ,NP3) ) 

F(5 ,NP2)= . 5*( . 5*(F(5 ,NP1)+F(5 ,NP2) )+F(5 ,NP3) ) 

F(6 ,NP2)=. 5*( . 5*(F(6 ,NP1)+F(6 ,NP2) )+F(6 ,NP3) ) 

810  UV(1)=-.5*ENU(1)*DUDDY(1)/RH0(1) 

VISCOS=VISMIX*DSQRT(TA)/RHOA 
IF (MODEL. EQ. 1)  GO  TO  812 
CE1=C0EFEP/DSQRT (RENO ) 

CE2=C0EFED/DSQRT(REN0) 

812  CE3=CE1 

DO  850  1=2, NP2 

UV(I)=-.5*(ENU(I)+ENU(I-1))*DUDDY(I)/RHO(I)*1.0 
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VV(I)=2.*(CC2-1.+CC1+(CM2/CC2-CC2)*BU0Y*UT(I)/FS(2,NP3) 

1  /F(5,I))*F(4,I)/(3 .*CC1) 

VT ( I ) = - 2 . *F (4,1) *VV ( I ) *DTDDY (I)/(F(5,I)*(2.*C1T 
1  +CT) ) 

IF(NALG.EQ.O)  GO  TO  820 

C - TT(I)  HERE  IS  ALG  -  TT  . 

TT(I)=-2.*F(4,I)*VT(I)*DTDDY(I)/(F(5,I)*CT) 

F(6,I)=TT(I) 

820  CONTINUE 

UT(I)=( -UV(I)*DTDDY(I) -VT(I )*DUDDY(I)*( 1 . +C2T)+BU0Y*F(6 , I ) 

1  *( 1 . -C2TM) /FS (2 ,NP3) )*F(4,I)/(. 5*F(5 , I)*(2*C1T 

2  +0 . 5*CT) ) 

PROK(I )=RHO(I)*( -UV(I )*DUDDY ( I)+BUOY*UT(I ) /FS (2 ,NP3) ) 
RDY=.5*R(I)*(Y(I+1)-Y(I-1)) 

IF  (I.EQ.2)  RDY=YI*R(2) 

IF  (I.EQ.NP2)  RDY=YE*R(NP2) 

BUPROK ( I ) =RHO ( I )*BUOY*UT ( I ) /FS ( 2 , NP3 ) 

SU(4 , I)=(PROK(I) -RH0(I)*1 . 00*F (5,1) )*RDY 
PROKM=PROK ( I ) 

IF  (LESSON. EQ. 2)  PROKM=RHO ( I ) * ( -UV ( I ) *DUDDY ( I ) +BUO Y*UT ( I ) 

1  /FS(2,NP3)*(CE3/CE1) ) 

FFF=DABS(F(5 , I) ) 

IF (MODEL. EQ.l)  GO  TO  830 

EPSPK=RDY*DSQRT (FFF/ VI SCOS ) * (CE l*PROKM) 

EPSDT=RDY*DSQRT(FFF/VISC0S)*(CE2*F(5,I)*RH0(I)) 

GO  TO  840 

830  EPSPK=RDY*FFF/F (4 , I )* (CE 1*PR0KM) 

EPSDT=RDY*FFF*(CE2*FFF) /F (4 , I )*RHO ( I ) 

840  SU (5 , I)=EPSPK-EPSDT 

SU (5 , I )=FACTOR*SU (5 , I ) 

PR0T(I)=RH0(I)*(-2.*VT(I)*DTDDY(I)) 

850  SU(6,I)=(PR0T(I)-RH0(I)*CT*F(5,I)*F(6,I)/F(4,I))*RDY 

F(4,2)=2.*(2.*F4KIN-F(4,1))-F(4,3) 
F(5,2)=2.*(2.*F5KIN-F(5,1))-F(5,3) 
F(6,2)=2.*(2.*F6KIN-F(6,1))-F(6,3) 
F(6,NP2)=2.*(2.*F(6,NP2)-F(6,NP3))-F(6,NP1) 
F(5,NP2)=2.*(2.*F(5,NP2)-F(5,NP3))-F(5,NP1) 
F(4,NP2)=2.*(2.*F(4,NP2)-F(4,NP3))-F(4,NP1) 

C . 

RETURN 

100  FORMAT (1H,A8 , 1P11E11.3/(9X,11E11.3)) 

END 

SUBROUTINE  STRIDE ( I SW) 


C/ . SUBROUTINE  FOR  PROGRAM  GENMIX  4A 

C/ .  D.B. SPALDING,  IMPERIAL  COLLEGE,  1972 


C/  THIS  SUBROUTINE  PERFORMS  THE  SAME  OPERATIONS  AS  THE  ONE  IN  GENMIX4A 
C  BUT  MORE  ECONOMICALLY.  THE  A,B,C  ARRAYS  ARE  ONE-DIMENSIONAL.  SOME 
C  OFTEN  USED  FUNCTIONS  OF  OM  ARE  STORED,  AND  A  D  ARRAY  SAVES 
C  UNNECESSARY  ARITHMETIC  IN  THE  TDMA  OPERATION. 

C . - . 

IMPLICIT 'DOUBLE  PRECISI0N(A-H,0-Z) 


on  o  no 
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DIMENSION  A2(7) ,ANP2(7) ,B2(7) ,BNP2(7) ,C2(7) ,CNP2(7) ,D2(7) ,DNP2(7) , 

1  AHLPT(83) ,B0MT3(S3) ,FDIFE(7) , 

2  FDIFI (7) ,GE(7) ,GI (7) ,PB0M(83) ,PGOM(83) ,THLPT(83) , TTFF ( 7 ) 

DIMENSION  B0M(83),0MP0M(83) 

COMMON/COMA/ A(83) ,AJE(7) ,AJI(7) ,B(83) ,C(83) ,CSALFA,D(83) ,DPDX(83) , 

1  DX,EMU(83) ,F(7,83),FS(5,83), I AX, IEND , IFIN, INDE (7 ) , INDI (7 ) , I  OUT, 

2  ISTEP,ITEST,IUTRAP,JS,JSW,JV, JY.KEX, KIN, KRAD,N,ND2,NF, NOVEL, NP1, 

3  NP2 ,NP3 ,0M(83) , OMD (83) ,P(83),PEI,PR(7), PREF(7 ,83) ,PSIE,PSII ,R(83) 
4,RHO(83) ,RME,RMI,RU(83),SD(7,83) ,SU(7,83) ,TAUE ,TAUI ,U(83) ,XD,XU, 
5Y(83) ,YE,YI ,ENU(83) ,NDEQ,BPI ,BPE ,DK1 (83) ,DK2(83) ,EDK1(83) ,EDK2(83) 
6 ,US(83) .FACTOR 

GO  TO  (1000,2000,3000,4000),  ISW 


1000  IF (I STEP)  1003,1003,1100 
1003  OMI=.5*OM(3) 

OME=.5*(l . -OM(NPl) ) 

DO  1002  1=2, NP2 
BOM(I)=OM(I+l) -OM(I-l) 
BOMT3(I)=3.*BOM(I) 

OMPOM ( I )=OM ( I )+OM ( 1+1 ) 
1002  OMD(I)=OM(I+l)-OM(I) 
OMD(l)=BOM(2) 

BPE=1 . 

BPI=1 . 


Y ( 1 ) =0 . 

IF(KRAD.EQ.l)  GO  TO  1100 
DO  1001  1=1, NP3 
1001  R(I)=1 . 

R25=l . 

RN15=1 . 

IF(ITEST.NE.O)  WRITE (6 , 9010)  (R(I) , 1=1 ,NP3) ,R25 ,RN15 
. - . .  CALCULATION  OF  RHO*U  'S  . 

1100  DO  1101  1=1, NP3 
IF(RH0(I) ,GT.0. )  GO  TO  1101 
WRITE (6 , 1108)  RHQ ( I ) , I , RHO ( 1 ) 

1108  FORMAT! 36H  **********  NEGATIVE  OR  ZERO  RHO (I )=, 1PE11 . 3 , 6H  AT  I=, 
1  13 , 6X, 21HSET  TO  ABS  OF  RHO( 1)=,E11 . 3 , 17H  ********  STRIDE 1) 
RHO(I)=DABS(RHO(l) ) 

1101  RU ( I ) =RHO ( I ) *U ( I ) 

RU3=RU(3) 

RUN1=RU(NP1) 

DO  1102  1=2, NP1 

1102  RU(I)=.5*(RU(I) +RU ( 1+ 1 ) ) 

IF(ITEST.NE.O)  WRITE(6 ,9010)  (RU(I) , 1=1 ,NP3) ,RUN1 ,RU3 ,PEI 

. .  CALCULATION  OF  Y  '  S  AND  R  '  S  . 

.  Y'S  FOR  PLANE  GEOMETRY 


YI=PEI*OMI/ (BPI*RU(2) ) 


Y(3)=YI+PEI*0M(3)/(RU(2)+RU3) 

Y(2)=2.*YI-Y(3) 

DO  1103  1=4, NP1 

1103  Y(I)=Y(I - l)+PEI*OMD (I-1)/RU(I-1) 

YN15=Y (NP1)+PEI *OMD (NP 1 ) / (RU (NP 1 ) +RUN 1 ) 

YE=PEI*OME/ (BPE*RU(NP1) ) 

Y(NP3)=YN15+YE 

Y(NP2)=2.*YN15-Y(NP1) 

IF(KRAD.EQ.O)  RETURN 

c  . -  Y'S  AND  R'S  FOR  AXI SYMMETRICAL  GEOMETRY 

IF (CSALFA . EQ . 0 . )  GO  TO  1110 

c  . . .  CSALFA  NE  ZERO 

COSD2=.5*CSALFA 

IF (R( 1) .NE . 0 . )  GO  TO  1105 

C .  R(1)=0. 

DO  1106  1=2, NP3 
Y(I)=DSQRT(DABS (Y(I) /COSD2) ) 

1106  R(I)=Y(I)*CSALFA 
YI=DSQRT(DABS (YI/COSD2) ) 

YN15=DSQRT(DABS (YN15/COSD2) ) 

GO  TO  1107 

c  .  R(l)  NE  0. 

1105  R1D2=.5*R(1) 

R1D2SQ=R1D2*R1D2 
DO  1104  1=2, NP3 

Y ( I ) =Y ( I ) / (R1D2+DSQRT (DABS (R1D2SQ+C0SD2*Y (I) ) ) ) 

1104  R(I)=R(1)+Y(I)*CSALFA 

YI=YI/ (R1D2+DSQRT(DABS (R1D2SQ+C0SD2*YI ) ) ) 

YN15=YN15/ (R1D2+DSQRT (DABS (R1D2SQ+C0SD2*YN15 ) ) ) 

1107  R25=R(1)+YI*CSALFA 
RN15=R(1)+YN15*CSALFA 
YE=Y (NP3) -YN15 
RETURN 

c  .  CSALFA  EQ  ZERO 

1110  DO  1111  1=2, NP3 
Y(I)=Y(I)/R(1) 

1111  R(I)=R(1) 

YI=YI/R(1) 

YN15=YN15/R(1) 

R25=R(1) 

RN15=R(1) 

YE=Y (NP3) -YN15 
RETURN ^ 

c  .  PRELIMINARIES  FOR  COEFFICIENTS 

2000  PX=PEI/DX 
PD8=.125*PX 
PD4=PD8+PD8 
G=RMI -RME 
ARMI=DABS(RMI) 

ARME=DABS(RME) 


o  o  o  n  o  o 


227 


GD4=.25*G 
PG=PX+G 
PGD8=. 125*PG 
PGD4=PGD8+PGD8 
RMID2=.5*RMI 
DO  2004  1=2, NP2 
PBOM( I ) =PX*BOM ( I ) 
2004  PGOM ( I ) =PGD4*0MD ( I ) 
P40MP=PD4*B0M(2) 


T1=0 . 


-  GRID  POINT  2 
TAUI ,  BPI,  T1 


IF(KRAD.EQ.O)  BPI=. 33333+ . 66667*RU( 1 )/RU(2) 

IF(KRAD.EQ.l)  BPI=(R(1)*(5 .*RU(1)+RU(2) )+3 .*R25* 

1  (RU( 1)+RU(2) ) )/6 . / (R(1)+R25)/RU(2) 

. - .  BOUNDARY  COEFFICIENTS  FOR  VELOCITY 

2002  HLP=RMID2-GD4*0MP0M(2) 

AHLP=DABS (HLP) 

THLP=HLP+HLP 
THLPT(2)=THLP 
TP=EM  U(2) 

TTP=TP+AHLP+DAB  S (TP - AHLP ) 

A(2)=TTP-THLP-T1-PG0M(2) 

B(2)=2.*T1+RMI+ARMI 

C(2)=P40MP*(3.*U(2)+U(3))-US(2) 

D(2)=A(2)+B(2)+PBOM(2) 

.  BOUNDARY  COEFFICIENTS  FOR  F'S 

IF(NF.EQ.O)  GO  TO  2304 
DO  2300  J=1 ,NF 

IF(J.GE.NDEQ+1 .AND. J.LT.4)  GO  TO  2300 
TPF2=TP/PREF (J, 2) 

TTPF ( J ) =TPF 2+AHLP+DAB  S (TPF 2 - AHLP ) 

T1F=0 . 

FDIFI ( J)=0 . 

2302  A2 ( J)=TTPF ( J) -THLP-T1F-PG0M(2)+. 5*SD( J, 2) 

B2 ( J)=2 . *T1F+RMI+ARMI 
SIMP=0 . 

IFCJ.EQ.5)  SIMP=SU(5,2)*(1. -FACTOR) /FACTOR 
D2 ( J)=A2 ( J)+B2 ( J)+PBOM(2) -2 . *SD ( J , 2) -SIMP 
T=-T1F*FDIFI ( J) 

TT=3.*F(J,2)+F(J,3) 

C2 ( J)=P40MP*TT+2 . * (T+SU ( J , 2 ) ) 

IF ( J . EQ . 5 ) C2 ( J) =P40MP*TT+2 . * (T+SU ( J , 2 ) ) 

2300  CONTINUE 


-  GRID  POINT  NP2 
TAUE ,  BPE ,  TNP3 

2304  TNP3=0 . 

IF (KRAD . EQ . 0 )  BPE= . 33333+ . 66667*RU (NP3 ) /RU (NP1 ) 
IF(KRAD.EQ. 1)  BPE=(R(NP3)*(5.*RU(NP3)+RU(NP1))+3.*RN15* 

1  (RU(NP3)+RU(NP1) ) )/6 . / (R(NP3)+  RN15)/RU(NP1) 

- -  BOUNDARY  COEFFICIENTS  FOR  VELOCITY 
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2310  HLM=RMID2-GD4*0MPQM(NP 1 ) 

AHL>1=DABS  (HLM) 

THLM=HLM+HLM 
TM=EM  U(NP1) 

TTM=TM+AHLM+DABS (TM-AHLM) 

P40MM=PD4*B0M(NP2) 

A ( NP 2 ) =2 . *TNP 3 - RME+ARME 
B (NP2 )=TTM+THLM-TNP3 -PG0M (NP 1 ) 

C (NP2 )=P40MM* ( 3 . *U (NP2 )+U (NP 1 ) ) -US (NP2) 
D(NP2)=A(NP2)+B(NP2)+PBOM(NP2) 

IF(NF.EQ.O)  RETURN 

C  .  BOUNDARY  COEFFICIENTS  FOR  F'S 

DO  2320  ,J=1 ,NF 

IF ( J . GE .NDEQ+1 . AND . J . LT. 4)  GO  TO  2320 
TMF=TM/PREF(J,NP1) 

TTMF=TMF +AHLM+DAB S (TMF-AHLM) 

TNP3F=0 . 

FDIFE ( J)=0 . 

2312  ANP2(J)=2 . *TNP3F- RME+ARME 

BNP2 ( J)=TTMF+THLM-TNP3F-PG0M(NP1)+ . 5*SD( J,NP2) 

SIMP=0. 

IF ( J . EQ . 5 ) . S IMP=SU ( 5 , NP2 ) * ( 1 . -FACTOR ) /FACTOR 
DNP2(J)=ANP2(J)+BNP2(J)+PB0M(NP2) -2 .*SD(J,NP2) -SIMP 
T=-TNP3F*FDIFE(J) 

TT=3.*F(J,NP2)+F(J,NP1) 

CNP2 ( J)=P40MM*TT+2 . * (T+SU ( J , NP2 ) ) 

IF(J.EQ.5)  CNP2(J)=P40MM*TT+2.*(T+SU(J,NP2)) 

2320  CONTINUE 
RETURN 

3000  DO  3005  I=3,NP1 
THLM=THLP 

HLP=RMID  2 - GD4*0MP0M ( I ) 

THLP=HLP+HLP 
THLPT(I)=THLP 
AHLP=DABS (HLP) 

AHLPT ( I ) =AHLP 
TTM=TTP 
TP=EM  U(I) 

TTP=TP+AHLP+DABS (TP-AHLP) 

A ( I ) =TTP -THLP - PGOM ( I ) 

B ( I ) =TTM+THLH - PGOM ( I - 1 ) 

C(I)=PD4*(B0MT3(I)*U(I)+0MD(I)*U(I+1)+0MD(I-1)*U(I-1))-US(I) 

D(I)=A(I)+B(I)+PBOM(I) 

D(I)=A(I)+B(I)+PBOM(I) 

3005  CONTINUE 

C .  . 

IF (KIN . EQ . 2 . AND . RU ( 1 ) .NE.O. )  U(1)=U(1)-DPDX(1)*DX/RU(1) 

IF (KEX . EQ . 2 . AND . RU (NP3 ) . ME . 0 . )  U(NP3 )=U (NP3 ) -DPDX(NP3 )*DX/RU (NP3) 
C .  SOLVE  FOR  DOWNSTREAM  U  'S  . 


B(2)=(B(2)*U(1)+C(2))/D(2) 


229 


A(2)=A(2)/D(2) 

DO  3048  1=3, NP2 
T=D (I)-B(I)*A(I-1) 

A(I)=A(I)/T 

3048  B(I)=(B(I)*B(I-1)+C(I))/T 
DO  3050  IDASH=2 ,NP2 
I=N+4-IDASH 
U(I)=A(I)*U(I+1)+B(I) 

3050  CONTINUE 

C  . 

IF (KIN. EQ. 3)  U(1)=.5*(U(2)+U(3) ) 
IF(KEX.EQ.3)U(NP3)=.5*(U(NP1)+U(NP2) ) 

C 

3013  IF(NF)  3060,3060,3014 

3014  DO  3320  J=1,NF 
IF(J.GE.NDEQ+1.AND. J.LT.4)  GO  TO  3320 

C-  — . . .  SOLVE  FOR  DOWNSTREAM  F  'S  . 

A(2)=A2( J) 

B(2)=B2(J) 

C (2)=C2 ( J) 

D(2)=D2(J) 

A(NP2)=ANP2(J) 

B(NP2)=BNP2(J) 

C(NP2)=CNP2(J) 

D(NP2)=DNP2(J) 

DO  3002  1=3, NP1 
TTMF=TTPF ( J) 

TPF=EM  U(I)/PREF(J,I) 

TTPF ( J ) =TPF+ AHLPT (I)+DABS(TPF-AHLPT(I)) 

A(I)=TTPF(J) -THLPT(I) -PGOM(I) 

B (I)=TTMF+THLPT(I-1) -PGOM(I-l) 

C(I)=PD4*(BOMT3 (I)*F ( J, I)+OMD(I)*F (J, I+l)+OMD(I-l)*F ( J, 1-1) )+ 

1  2.*SU(J, I) 

SIMP=0 . 

IF ( J . EQ . 5 )  S IMP=SU ( 5 , I ) * ( 1 . - FACTOR ) /FACTOR 
3002  D(I)=A(I)+B (I)+PBOM(I) -2 .*SD(J, I) -SIMP 


B(2)=(B(2)*F(J,1)+C(2))/D(2) 

A(2)=A(2)/D(2) 

DO  3148  1=3, NP2 
T=D(I) -B(I)*A(I-1) 

A(I)=A(I)/T 

3148  B(I)=(B(I)*B(I-1)+C(I))/T 
DO  3150  IDASH=2 ,NP2 
I=N+4-IDASH 

3150  F(J,I)=A(I)*F(J,I+1)+B(I) 

C- .  ADJUST  F(J, 1)  AND  F(J,NP3) 

GO  TO  (3220,3220,3230) , KIN 
3230  F(J,1)=.5*(F(J,2)+F(J,3)) 

3220  GO  TO  (3320,3320,3330) ,KEX 
3330  F ( J ,NP3)=. 5*(F ( J,NP1)+F( J ,NP2) ) 
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3320  CONTINUE 


3060  XU=XD 


ITEST=0 

PSII=PSII-RMI*DX 
PS IE=PS IE -RME*DX 
PEI=PSIE-PSII 
ISTEP=ISTEP+1 
RETURN 


-*  STRIDE4 

4000  CONTINUE 
ND2=N/2 
NP1=N+1 
NP2=N+2 
NP3=N+3 
OM ( 1 ) =0 . 

0M(2)=0. 

0M(NP3)=1. 

ISTEP=0 

IEND=10000 

IAX=10000 

IOUT=10000 

DX=1 .E-30 

IFIN=0 

DO  4001  J=l, 7 
DO  4001  1=1, NP3 
SU(J,I)=0. 

4001  SD(J, I)=0 . 

RETURN 

9010  FORMAT ( 1H  ,1P11E11.3) 

END 

SUBROUTINE  PLOTS  (X, IDIM, IMAX.XAXIS , Y, JDIM, JMAX, YAXES , SYMBOL) 
C*****14.8, 

C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


SUBROUTINE  FOR  PLOTTING  J  CURVES  OF  Y(J,I)  AGAINST  X(I) . 

X  AND  Y  ARE  ASSUMED  TO  BE  IN  ANY  RANGE  EXCEPT  THAT  NEGATIVE  VALUES 
ARE  PLOTTED  AS  ZERO. 

X  AND  Y  ARE  SCALED  TO  THE  RANGE  0.  TO  1  BY  DIVISION  BY  THE  MAXIMA, 
WHICH  ARE  ALSO  PRINTED. 

IDIM  IS  THE  VARIABLE  DIMENSION  FOR  X. 

IMAX  IS  THE  NUMBER  OF  X  VALUES. 

XAXIS  STORES  THE  NAME  OF  THE  X-AXIS. 

JDIM  IS  THE  VARIABLE  DIMENSION  FOR  Y. 

JMAX  IS  THE  NUMBER  OF  CURVES  TO  BE  PLOTTED,  (UP  TO  10). 

THE  ARRAY  YAXES (J)  STORES  THE  NAMES  OF  THE  CURVES. 

THE  ARRAY  SYMBOL(J)  STORES  THE  SINGLE  CHARACTERS  USED  FOR  PLOTTING. 


* 


* 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  X(IDIM) ,Y( JDIM, IDIM) , YAXES (JDIM) , SYMBOL (JDIM) , 
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1  A (101) ,YMAX(20) 

DATA  DOT , CROS  S , B  LANK/ 1H . , 1H+ , 1H  / 

C*****  SCALING  X  ARRAY  TO  THE  RANGE  0  TO  50 
X>L4X=1 .  E-30 
DO  1  1=1 , IMAX 

1  IF(X(I) .GT.XMAX)  XMAX=X(I) 

DO  2  1=1, IMAX 
X(I)=X(I)/XMAX*50. 

2  IF(X(I ) . LT. 0 . )  X(I)=0. 

C*****  SCALING  Y  ARRAY  TO  THE  RANGE  0  TO  100 
DO  3  J=1 , JMAX 
YMAX( J)=l .E-30 
DO  4  1=1, IMAX 

4  IF(Y(J,I) .GT.YMAX(J))  YMAX(J)=Y(J, I) 

DO  3  1=1, IMAX 
Y(J,I)=Y(J,I)/YMAX(J)*100. 

3  IF(Y( J , I) . LT. 0 . )  Y(J, I)=0 . 

C*****  IDENTIFYING  THE  VARIOUS  CURVES  TO  BE  PLOTTED 
WRITE (6, 103)  XAXIS.XMAX 
WRITE (6 , 100)  (YAXES(I) , 1=1, JMAX) 

WRITE (6 , 106)  (SYMBOL(I) , 1=1 , JMAX) 

WRITE(6, 102)  (YMAX(I) , 1=1 , JMAX) 

DO  5  1=1,11 

5  A(I)=0.1*FL0AT(I-1) 

WRITE (6 , 101)  ( A ( I ) , 1=1 ,11) 

C*****  MAIN  LOOP.  EACH  PASS  PRODUCES  AN  X-CONSTANT  LINE 
DO  40  1=1,51 

IF ( I . EQ . 1 . OR . I . EQ . 5 1 )  GO  TO  32 
GO  TO  33 

C*****  ALLOCATE  .  OR  +  AS  MARKER  ON  THE  Y-AXIS 

32  DO  30  K=l, 101 

30  A(K)=DOT 

DO  31  K=ll , 101,10 

31  A(K)=CROSS 

C *****  ALLOCATE  .  OR  +  MARK  ON  THE  X-AXIS,  ALSO  THE  APPROPRIATE  X  VALUE 

33  A(l)=DOT 
A(101)=DOT 
K=I-1 

46  K=K-5 
IF(K)48,47,46 

47  A(l)=CROSS 
A(101)=CROSS 

48  XL=0 . 0 2*FLOAT ( I - 1 ) 

C-w.***  CHECK  IF  ANY  Y(  X(I)  )  VALUE  LIES  ON  THIS  X-CONSTANT  LINE 
C*****  if  YES  GO  TO  41,  OTHERWISE  GO  TO  42 
DO  43  K=1 , IMAX 

IF(IFIX(X(K)+1 . 5) -I )  43,41,43 
C*****  LOCATE  Y (  X ( I )  ) 

41  DO  44  J=l, JMAX 
NY=Y ( J , K ) + 1 . 5 
A(NY)=SYMBOL(J) 
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44  CONTINUE 
GO 'TO  42 
43  CONTINUE 

C *****  PRINT  X-CONSTANT  LINE 

42  WRITE (6 , 105 )  XL, (A(K) ,K=1 , 101) ,XL 
C*****  PUTTING  BLANKS  INTO  X-CONSTANT  LINE 
DO  49  K=1 , 101 

49  A(K)=BLANK 
40  CONTINUE 

DO  50  1=1,11 

50  A ( I ) = . 1*FL0AT ( I - 1 ) 

WRITE(6 , 104)  (A(I) , 1=1 ,11) 

RETURN 

100  FORMAT (11H  Y-AXES  ARE ,5X, 10(1X,A10) ) 

101  FORMAT ( 1H0 , 2X, 11F10 . 1) 

102  FORMAT (15H  MAXIMUM  VALUES , 1P10E11 . 3) 

103  FORMAT (11H0X- AXIS  IS  ,A8,17H  .MAXIMUM  VALUE  =,1PE10.3) 

104  FORMAT (3X, 11F10. 1/1H1) 

105  FORMAT (2H  X,F6 . 2 ,3X, 101A1.F6 . 2) 

106  FORMAT ( 7H  SYMBOL, 11X, 10 ( IX ,A10) ) 

END 
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